es 
ou nc 


sy 
7 J the h “tives 
Aa fy ay) Basha eh Sieh, oo i i aS , i Rate Peace Fen he 
Way feta AUN A ta ALAS. Baek Ts hs Tine 
; : Aare eee eres EARS ANON TIA A i ween % a tat. ations “RMI LAs 
. J ‘ ; i 4 " 5 " tine. ne nA 5 fy SS KA te te Rae As 2 eae 7 ts Os Weta Sri eek ety a 
j : : F " u x p : ay gy ea he 3 rn esse SAR kree ee RRS 
F fi . of ‘ . ' Ae ‘ } iy Tae Te St Ke He, MRT TS BAB ac My Ta a Ps FTE, 
' ‘ : ; i 4 i eR Hatta i cit ; 4 te cle tel Ba ag hie Ib eas 


pe ae We Fy Ny NH 
ty Ahn hots Vas WT: Me sch te & ® Fate Nec shee ys AAA we Teresi 
any Tecte'nts Maita Sui ate Ae ae ea aire Se ak 
fe Ae lke ae ae Tas efit ate Tea NS 
HS be: E hae 42 AA 
2 . as a1 Tatas Tn ta The TH 
ak “ Ye VATS AS EIR Fe Be RT veh. hes ek ei att 
a : 4 A Bit AO ym ON ORE. Wie et Reeser tse etic Ve 
sg : d ' ‘ s he *, ne i i ax ta Meath Harte ¥ Me the ah Ree 
ns . i i ie PN REGS At ihe A Snetricer ys 
Hoke = Fes Pal Retest te 
be hots potrkk et Shem tate ND 
es Mee SEAT Tina Fn A Mote h Be Mace he Mt i 
Ne tee Fe Peete ae a thee Se AAS ta ah a Ae uc 
r Reh Res Ea ae % 
‘s % AS Th ACL. %. 
Re DN Tes He TT * 7 Be BY he Be he 
yw fag WN fhe BE Soe te 
et Witte: a Be 


yin ei 


Tek Aye Be, 1 a rea 
vis iohe Ry. SOW 


: . i fy ns hy hs. a earg Oe “ * shi ct Tas icy athe Ths 
, y sii: hs ae 4 “ ok Te, ‘ “Hes hi x : ra . te mes Ft A Aa 
$ : “ : A nN Z i as : A hse ee MT 
, a ch ec Rc 
Fis Ye eye Fos 
ah fe Ra RR 


For Reference 


NOT TO BE TAKEN FROM THIS ROOM 


“us 


Gb] eed Paw 


PCI O99 AF 9 oF 
Cee atu ee ae aad 


in das rested. pede ie Sel 
reed Prien oars 
ree aa s 





Pe Ruessen ie ein 
3 : ‘i i 3 $i f Wie Egeret As if Bee 
ES By, ee BF L f 4 fe Te ot aa PPI Aer pane i 


- AIBRIS 
NIOKRSUTATIS 
ne RUAEASTS 





Digitized by the Internet Archive 
in 2019 with funding from 
University of Alberta Libraries 


_https://archive.org/details/Kasper1981 








To beON LViBE Seen Obes ep caen 
RELEASE FORM 


NAME OF AUTHOR Bryon Lynn Kasper 
geo bee Pele Shs An Interference-Rejecting Radiometer for 
Low Frequency Astronomy 
DEGREES EORSWHICH THESIS WAS PRESENTED “Doctor of "Philosophy 
Teak omDEGREE GRANTED January 1981 
Permission is hereby granted to THE UNIVERSITY OF 
ALBERTA LIBRARY to reproduce single copies of this 
thesis and to lend or sell such copies for private, 
scholarly or scientific research purposes only. 
The author reserves other publication rights, and 
neither the thesis nor extensive extracts from it may 
be printed or otherwise reproduced without the author’s 


written permission. m 





em Vien Sue (Ore ae Brn A 


An Interference-Rejecting Radiometer for Low Frequency 
Astronomy 


by 





Bryon Lynn Kasper 


Petite SS 
Soe eee ee OP he PACUL TY UF “GRADUATE STUDDES AND RESEARCH 
Hee rh en PUP PEMENT OF THE REGULREMENTS FOR THE DEGREE 


OF Doctor of Phi losophy 


Electrical Engineering 


EDMONTON, ALBERTA 


SPRING 1981 





we 
« — £ 2! a *| ; - 5 = F] ex? 4 » 
Ai WO 97 1sssmofhsa Onksaslsh ssnete tae Fr. 


noriotrhdaas 


(ele WINMLNMS eS Taye fdbiaie renege 
PACUMi@OnmGRADUATESTUDIESSAND RESEARCH 


The undersigned certify that they have read, and 
recommend to the Faculty of Graduate Studies and Research, 
for acceptance, a thesis entitled An Interference-Rejecting 
Radiometer for Low Frequency Astronomy submitted by Bryon 
Lynn Kasper in partial fulfilment of the requirements for 


the degree of Doctor of Philosophy. 





Dedicated to my wife Vicky 


and to my Parents 


For years of encouragement and understanding. 





yer Fy 2 
| 
i= 


.gnibagraisbeu bes. fi rice 1 ey, 
a reabapucene, Yo sey 109 | 
t . | | 


i 





Abstract 
A serious difficulty for low frequency radio astronomy 
either on the earth’s surface or in any unshielded location 
near the earth is the presence of interfering terrestrial 
radio transmissions. Generally, however, interference will 
not exist simultaneously at all frequencies. Gaps in the 
spectrum which do not contain significant interference may 
be used to allow astronomical observations at times when 
adjacent interference would normally make observing 
impossible. 

An experimental system is designed and built to 
investigate the possibility of on-line detection and 
rejection of narrowband terrestrial interference. Digital 
cross spectral analysis via the fast Fourier transform 
allows narrowband interference to be identified as peaks in 
the spectrum. A low-cost FFT processor designed and built 
for this purpose is described. The processor is capable of 
calculating a 256-point FFT in 2.458 msec for a real-time 
bandwidth of 52 kHz. 

The expected distributions of cross and auto spectral 
components in the absence of interference are derived. 
Robust estimation which weights outlying points less heavily 
than those near the center of a distribution is employed to 
reject interference and estimate the centers of the cross 
spectra. The first deciles of the auto spectra are used as 
simple but accurate estimates of scale for the cross 


spectra. A robust procedure combining outlier rejection, 
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Huber M-estimation, and biweight M-estimation then 
determines the correlated broadband noise levels of the 
cross spectra while minimizing the effects of outliers. 
Experimental testing of the system using an 
interferometer at 22.25 MHz during the winter of a solar 
maximum is described. Interference removal is found to be 
highly successful during periods with less than 30% and 
occasionally up to 50% of the bandwidth being occupied by 
interference. Low level interference during the night was 
encountered on 13 of 15 nights and was eliminated 
completely. Additional observing time of from 60 to 90 
minutes was generally obtained in both the mornings and 
evenings. The distribution of interference amplitudes was 


found to be reasonably well represented by a power law. 
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1. Introduction 
Nearly all of man’s Knowledge of the universe beyond our own 
planet has been derived from the study of electromagnetic 
radiation emitted by the cosmos. Astronomy began at visible 
wavelengths and was confined there until 1932 when Jansky 
discovered extraterrestrial radio emissions and opened the 
door to an exciting new field of radio astronomy. Presently, 
astronomy is conducted from radio wavelengths up to gamma 
rays, with additional studies of cosmic rays, gravitational 
waves, and neutrino emissions. 

Radio astronomy had its beginnings at decametre 
wavelengths, but quickly moved to higher frequencies as 
radio technology progressed. Because of a number of 
difficulties mainly associated with the ionosphere, low 
frequency astronomy has received comparatively little 
attention. In recent times, new regions of the spectrum have 
often yielded new and and unexpected JHimaower Ie. A thorough 
exploration of the sky at decametric and longer wavelengths 
could therefore add considerably to our understanding of the 
universe. 

The subject of this thesis is an experimental attempt 
to overcome one of the problems which plagues low frequency 
astronomy, namely interference from terrestrial radio 
transmissions. Digital spectral analysis is employed to 
produce spectra in which narrowband man-made signals stand 
out from broadband cosmic signals. By applying robust 


estimation techniques to the spectra, it is possible to 
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exclude the narrowband, interfering signals and accurately 
measure cosmic signals. 

It is hoped that through the removal of terrestrial 
interference it will be possible to considerably enhance the 
accuracy Of observations and to greatly extend the amount of 
time during which suitable observing conditions are present 
for low frequency astronomy. 

This chapter will discuss decametric astronomy and the 
difficulties which it presents. The ionosphere and its 
effects on terrestrial and astronomical radio signals will 
be described. Previous attempts to deal with the problem of 


terrestrial interference will also be discussed. 


1.1 Decametric Astronomy 

The longest wavelengths where astronomy has been 
attempted to any degree are in the decametric region. There 
are a number of areas of astronomy for which decametric 
observations are important and can contribute significantly 
to the understanding of astrophysical phenomena. The 
following discussion is based upon Routledge [1], Dewdney 
U2wand Kraus 1/3); 

Decametric astronomy may be divided into studies of 
discrete sources and extended sources. Discrete sources are 
point-like sources which are smaller in angular extent than 
the resolving power of a telescope, for example some distant 
Galaxtes. intensittes of discrete sources are expressed in 


PeGmsmote total Coserved flux density, S(T}, with the unit 
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being the Jansky (10-28 watt m-2 Hz-1'). Extended sources are 
larger in angular extent with an intensity which changes 
smoothly over the source, for example large interstellar 
clouds of ionized hydrogen (HII). Measurements are in terms 
of brightness, B(f), in Jansky-steradian-', or equivalent 
black-body temperature, T(f), in degrees Kelvin. 

One of the most important features of radio sources is 
the spectrum (intensity vs frequency). No strong spectral 
lines have been found below 1000 MHz, thus at low frequencies 
radio emissions are not frequency selective. Instead, al] 
frequencies are present producing a continuum of radiation 
with no discontinuities. Intensity is found to vary with 
frequency in a smooth fashion, 


S(£) a(f) 


i} 


(const) £. 


(const) 784) 


MW 


or T(£) 


a(f) and 8(f) are Known as the flux density spectral index 
and temperature spectral index, respectively. They are 
generally positive, implying that intensity increases as 


frequency decreases. 


1.1.1 Discrete Sources 

Discrete sources may be classified as galactic or 
extragalactic. Galactic sources include supernova remnants 
and small HII regions. 

Supernova remnants radiate due to the synchrotron 


mechanism. A supernova explosion produces an expanding cloud 
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of matter containing relativistic electrons (about 1012 eV) 
and magnetic fields. The electrons, which may continue to be 
generated by remnants of the explosion, interact with the 
magnetic field to produce synchrotron emission. If the 


electrons have an energy spectrum of the form 


N(E) = (const) ER 


where N(E) = number of electrons as a function of energy 
E = energy of electron 
v = energy spectrum index 

then the synchrotron radiation will have a characteristic 


spectrum given by 


Gp ow Oe 


Observations have shown that the above spectrum is 
typical of many non-thermal radio sources over a wide range 
of frequencies. However, at low frequencies the intensity 
would tend towards infinity, thus at some point a turn-over 
in the spectrum is necessary. Low frequency astronomy is 
necessary for investigations of the turn over of non-thermal] 
radio sources, in order that the nature of such sources may 
be more clearly understood. 

The study of synchrotron emissions produces information 
about the physical structure and the evolution of the 
sources. Calculations may be made of plasma densities, 
magnetic field strengths and relativistic particle 


densities. A number of mechanisms believed to cause 
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deviations from a straight power-law spectrum can also be 
examined. These include variations from a power law in the 
source of relativistic electrons, multiple sources of 
electrons with different energy spectra, thermal absorption 
by ionized gas either within the source or between the 
source and the earth, and synchrotron self-absorption within 
the source. 

An interesting class of sources which has been 
discovered through decametric astronomy is steep spectrum 
sources possessing higher spectral indices than expected at 
low frequencies. Study of these objects, for example the 
source at the center of the Crab Nebula, may produce clues 
to the history of supernova remnants. 

Extragalactic sources are generally galaxies or 
clusters of galaxies, some of which emit enormous amounts of 
energy. These extremely powerful radiation sources have 
attracted much attention because they defy many previously 
accepted notions about the universe. Optical counterparts of 
strong radio sources include giant elliptical galaxies, 
N-type galaxies, and quasi-stellar objects (quasars). 

Extragalactic sources are non-thermal in nature and 
often emit because of the synchrotron mechanism. Low 
frequency studies of such objects may help to explain the 


origin of the enormous energies produced. 


1.1.2 Extended Sources 


With the comparatively low resolution of present 
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decametric telescopes, observable extended sources are 
confined to our own galaxy. The major type of extended 
source is interstellar plasma produced by low-intensity 
ultraviolet radiation from early-type stars and by 
low-energy cosmic rays. Although ionized hydrogen does emit 
thermal radiation at decametric wavelengths, it is most 
easily observable because of its absorption of non-thermal 
emissions from the galactic background. Free-free absorption 
causes HII regions to appear as localized areas of reduced 
brightness or reduced spectral index against the brightness 
of the galaxy. Low frequency telescopes provide one of the 
most sensitive methods of detecting interstellar HI1. 
Decametric mapping of the distribution of HII in the 
galaxy would be an important contribution to studies of 


galactic structure and galactic dynamics. 


1.1.3 Existing Decametric Telescopes 

High-quality low frequency radio telescopes are 
difficult and expensive to build. The frequencies and 
angular resolutions of 21 telescopes constructed since 1958 
for use below 100 MHz are summarized in a CCIR Report [4]. 
Below 20 MHz the best resolution obtained has been a few 
degrees, and major studies have been limited to two source 
surveys ([5], [6]) and measurements of the galactic 


background ([7], [8]). 
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1.2 Limitations to Low Frequency Astronomy 
There are three major factors which tend to discourage 
attempts at low frequency astronomy: 

1. Very large antenna dimensions (several Kilometers or 
more) are required for good resolution. 

2. Low frequency signals are highly susceptible to 
distortion by the ionosphere. The ionosphere can cause 
refraction, scintillation and absorption of signals from 
cosmic sources. 

3. No frequencies below 20 MHz have been set aside 
exclusively for radio astronomy'!. The radio spectrum in 
this region is used extensively for terrestrial 
communications. Long distance propagation via ionospheric 
reflections is possible, making it extremely difficult if 
not impossible to find sites for telescopes which are 
immune to terrestrial interference. 

The resolution of an antenna in a plane is about equa! 
to the half-power beam width in that plane [3]. In an 
idealized situation, the number of resolvable sources 


distributed uniformly over the sky is approximately given by 


where &, = antenna beam solid angle, rad*. Beam solid angle 


is related to antenna effective aperture, A., and wavelength 


weer www eX’ ek— MX MX KH WX MX— MK MY VM MX MX 


1A new band from 13360 to 13410 kHz was allocated to radio 
astronomy during the 1979 World Administrative Radio 
Conference in Geneva. This allocation is not exclusive but 
is to be shared with the Fixed service. 
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as follows: 


= 2 
n, d [A 


To maintain a given resolution the linear dimensions of an 
antenna must increase proportionally with wavelength. For 
example, to obtain a resolution of 5 arc minutes, as is 
possible at microwave frequencies with large dish antennas, 
would require an antenna 10 km in extent at 20 MHz. No 
antennas of that size presently exist. The largest ones 
built to date are about 3 km. 

Fortunately, an antenna 10 km in extent or 100 km?2 in 
area need not be completely filled to produce the desired 
resolution. Early telescopes (for example Mills [9]) were 
built in the form of crosses, T’s and circles to take 
advantage of this fact. Development of the principles of 
aperture synthesis by Ryle and Hewish [10] showed that an 
antenna may be decomposed into a number of elemental units, 
and that it is only the relative positions or spacings of 
the units which are fundamentally important. Measurements 
with a resolution equivalent to that of a very large filled 
antenna may be made by combining signals from smaller 
antennas with all the elemental spacings between them. 

For sources which are not time-variable, the 
observations at different spacings need not be made 
simultaneously. Therefore an antenna of any chosen size wey 
be synthesized by using just two small elementary units 


which can be moved to all of the required positions to 
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produce the spacings of the desired large antenna. Indeed, 
it is possible to let the rotation of the earth provide part 
of the movement of the two antennas relative to the sky. 
Such earth rotation aperture synthesis as described in 
Fomalont [11] is now commonly used for high resolution 
mapping of the sky. 

The result of all this is that resolution at low 
frequencies is not fundamentally limited by the large 
antenna sizes required. Limitations due to the ionosphere 


are much more difficult to overcome. 


1.3 The Ionosphere 

The following discussion is derived mainly from Davies 
Gil 

The earth’s ionosphere is produced by high-energy 
radiation from the sun, mainly in the ultraviolet and soft 
X-ray regions. Electron density depends upon many factors, 
including the intensity and spectrum of incoming radiation 
(some of which may be screened by higher ionospheric 
layers), atmospheric density, and chemical composition at 
the altitude in question. Solar control results in a strong 
dependence on the time of day (or night), the season, 
latitude, and the level of solar activity at a given time. 
Also, the electricel nature of the ionosphere causes it to 
interact a great deal with the earth’s magnetic field. 

The ionosphere has attracted a tremendous amount of 


study since the 1940's, largely due to its importance in 
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radio communications. Although most of its properties have 
been characterized and its long-term behavior is reasonably 
well understood, day-to-day ionospheric activity is highly 
variable and somewhat unpredictable, much like the earth's 
weather. 

The ionosphere is ne ner etd divided into a number of 
regions, the most important of which are termed the D, E and 
F regions. The D region extends in height from approximately 
50 km to 90 km, the E region from 30 km to 130 km, the F 
region from 130 km to perhaps 500 km. 

D-region ionization is largely a daytime phenomenon, as 
the relatively high atmospheric density allows rapid 
recombination after the ionizing radiation disappears. The D 
region is mainly responsible for the absorption of radio 
signals. AM radio transmissions, for example, are limited in 
daytime propagation distance by D-region absorption, whereas 
at night they can propagate much farther. 

The E region is also normally significant only in 
daylight. It reflects radio signals below a maximum of about 
Re ZeerneeUSUdl lah eCTOMelOnmiZzali1on VS not important for 
terrestrial propagation above this frequency. However, there 
is a version of the E region Known as sporadic E— which is 
sometimes very significant. Sporadic E, as the name 
Slidge sts maulcudaliaicicm nmuUnpredi1clable and Often sloca lized 
phenomenon. When it is present, sporadic E can allow much 
longer distance radio communication at higher frequencies 


than would normally be expected. 
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Durinceetnencay the Fureqion splits into two ayers. The 
lower one is termed the 4 layer and the upper one, the F 
layer. The Fy layer, which disappears at night, reflects 
radio frequencies somewhat higher than those reflected by 
the E region. Like the E region, the electron content of the 
Fy layer is higher in summer than in winter. 

The part of the ionosphere which is the most important 
for HF radio communication, and also the most detrimental to 
radio astronomy, is the F. layer. This layer almost always 
contains the maximum electron density in the ionosphere. 
Recombination occurs very slowly, hence the Fo layer is 
present at night and reaches a minimum electron density just 
before sunrise. The seasonal minimum is during the winter 
night, as it is for other layers. Surprisingly, however, the 
seasonal maximum occurs during the winter day rather than 
during the summer day. This peculiarity of the aS layer is 


termed the "winter anomaly’". 


1.4 Ionospheric Radio Propagation 

The refractive index of an ionized medium depends upon 
the electron density. If the collision frequency between 
electrons and neutral molecules or ions is low (as in the 
rarified upper atmosphere of the E and F regions) then 
absorption of radio energy is negligible. If the magnetic 
field is assumed to be zero, then the refractive index is 
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where fae is the plasma frequency given by 


Ne? bisa 
fy = ae (Elpe2s) 


with N = electron density 


electronic charge 


e 
Ena ts permittivity of free space 

m = electron mass 

If a magnetic field is present the situation becomes 
more complex. Two electromagnetic waves which are circularly 
polarized in opposite directions (termed “ordinary” and 
"extraordinary" waves) will propagate with two different 
indices of refraction. The refractive index for the 
extraordinary wave is the higher of the two. 

The frequency at which the refractive index becomes 
zero is termed the critical frequency, fe: At this frequency 
or below it a normally incident wave cannot propagate but is 
instead entirely reflected. Extraordinary waves exhibit a 
higher critical frequency, Fen} than do ordinary waves, fe, 
due to the difference in refractive indices. 

At frequencies above i an obliquely incident wave may 
still be completely reflected. The maximum frequency which 


will be completely reflected at a given angle of incidence 


o. is given by 


A 
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In propagation studies the term maximum usable frequency 
(MUF) is used for the maximum frequency which can be 
reflected between two points a given distance apart on the 
earth's surface. The angle of incidence for reflection at a 
given distance depends upon the height of the reflecting 
layer. If the curvature of the earth and ionosphere plus the 
normal height of the Fs layer are considered, it may be 
shown [18] that the largest angle of incidence possible for 
reflection back to the earth’s surface is usually about 74 
degrees. The maximum usable frequency at this limiting angle 
is 


MUF = f sec 74° = 3.6 f esa) 


From a propagation viewpoint, maximum usable frequency 
is quite important. D-region absorption increases with 
wavelength, so it is advantageous to operate on the highest 
frequency which will support reliable reflections. HF radio 
users, with the help of ionospheric propagation forecasts, 
often change operating frequencies to obtain maximum 
propagation distances. 

Critical frequency and hence MUF vary regularly through 
the day as the sun’s angle changes and through the seasons 
as the length of day and solar angle vary. In addition, a 


strong dependence exists on the 11-year sunspot cycle. 


13 
















ball rT on PG 
> a 
2 = WY 
An 3 = 4 ry 
7 ar | 
es Ne ‘on Med g3 BetOuls no? Agee 


a) riety 18 1? Kinks, anz 107 beau 





> foes fsreib revem s @e7ateq.owl neewsed be 
? y 1 - 
Po 
¢ ; —, | ie el) “4 a 
Oi 75° 3 10 sSNeoi an ‘o Sipns<or 2087 we ent: 
ifositsn sol tc pian eT nequ Sepregen soraiath 
t « ‘7 %o 23'it) TLS <)| 
T 18 Velo, Sm to 2fipt 
z 4 — ; 
‘ . > a = | 
i eB | tt) soi ro? 
j age él dseunnbnisaey ear 
‘ 7 
‘ se > GR -) 
= if; GT , af Wwery. Mcp sa4 2HGO1G; 5 
Alsw. 2 g TOT hay se narpse >=) 5 7a [ott 
A earl? 9 ShETOq iT Aaurepets BvES 270 39 “oe, a orrat« 
ie. i Rates 
ec 7 Py a r ° pe ‘ hao & * o. 
: y? +8. ss aa | ways Ped Lf WwW fri, id -.) 


fh227? nOorieqeno7 3 SSrQaorsy to ; if sort ast? Ante tl 

nunca ces do. oy aroneraes 1 co saan ara 
| i” ‘ Avie ve, -esonadatb ne C 
>) ange! bis, i C eat a 


a2cae art riquio sf — gS apMSM2 signs 
at vi Ps —: 4 


. nai Son's ete shone as 


on 





A? vi se byes cna wile 


Figure 1.1, taken from Jordan and Balmain [19], illustrates 
diurnal critical frequency variations for the E, Pas and F, 
layers in winter and summer, and at sunspot minimum and 
maximum. The critical frequency of the F, layer, tee oe is 
higher during the day than at night, and higher in winter 


than in summer. At sunspot maximum LP a is approximately 


twice its value at sunspot minimum. 


1.5 lonospheric Effects on Radio Astronomy 
The major effects which the ionosphere has on 
astronomical signals are absorption, refraction and 


Sci aisiia uO. 


leo leeAOSOrD tion 

Virtually all absorption at decametric wavelengths 
occurs in the D region which is present only during the day. 
The degree of absorption is dependent on frequency and 
reaches a broad maximum near the gyromagnetic frequency for 
an electron in the earth’s magnetic field (1.2 to 1.4 MHz). 
A device called a riometer [20] may be used to measure 
absorption by monitoring the average signal strength from a 


broad area of the sky. 


1.5.2 Refraction 
For a plane, stratified ionosphere, incoming waves at a 
frequency f above the F-layer critical frequency Fee will 


propagate through the ionosphere provided that the source is 
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Figure 17-8. Diurnal variation of critical frequency and virtual height of 
the regular ionospheric layers: (a2) summer at period of sunspot minimum; 
(b) summer at sunspot maximum; (c) winter at sunspot minimum; (@) 
winter at sunspot maximum. 


Figure 1.1. Diurnal Critical Frequency Variations 
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within a cone centered on the zenith and defined by a 


limiting zenith angle [4] 


Ds = Giee Cos ) 


For a spherical ionosphere the actual limiting angle will be 
slightly larger. Spherical stratification produces a slight 
refraction of signals passing through the ionosphere such 
that the apparent zenith angle is less than the true zenith 
angle by a small amount [21]. 

Uniform horizontal variations in -ionospheric electron 
density, such as occur at sunrise and sunset, are another 
cause of refraction [22]. Such gradients may cause 
refraction of a few degrees if the observing frequency is 
less than twice the critical frequency. 

Correction for refraction may be accomplished by 
monitoring the apparent positions of calibration sources, 
provided one is using an antenna which can provide a narrow 
beam to track the sources at any given time. However, as 
refraction is a time-varying phenomenon, it would be 
extremely difficult to correct for in the case of aperture 


synthesis which relies on time-invariant source positions. 


1.5.3 Scintillation 

A transient phenomenon which can seriously distort 
signals traversing the ionosphere is scintillation. If the 
ionosphere contains irregularities or inhomogeneities in 


electron content, then the index of refraction will vary for 
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rays having different paths. As a result, a plane wave 
impinging on the ionosphere will have horizontal phase 
variations when it emerges below the ionosphere. Such phase 
variations are converted to both phase and amplitude 
variations due to interference effects as the wave 
propagates towards the earth’s surface. As a result, a radio 
source appears to "twinkle" just as a star does optically 
when seen through the turbulent atmosphere. A discussion of 
scintillation is contained in Briggs [23]. 

The irregularities which cause scintillation occur 
mainly in the F region and have been correlated with 
observations Known as "spread-F" during ionospheric sounding 
[14]. Spread-F consists of multiple reflections of probing 
radio signals from points within the F region, and is 
consistent with the idea of irregularities in electron 
density in this region. 

Fortunately, scintillation is present only on occasion. 
It is most frequently observed near the equinoxes of a solar 
maximum [15]. During solar minimum, scintillation occurs 
mainly at night, whereas during solar maximum it is equally 
likely by day or night. Geographically, scintillation is 
least likely at temperate latitudes. For astronomy, the 
severity of scintillation is dependent upon the scale of 
irregularities and upon antenna size. Irregularities are 
generally found to be from about 0.5 to 10 km in extent [4]. 
For antennas larger than 2 km, amplitude and phase 


variations will be produced across the aperture and may 
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distort the antenna beam. For smaller antennas (and larger 
beam widths) larger phase variations can be tolerated. 

The effects of scintillation are most severe at low 
frequencies because the ionosphere’s refractive index goes 
from 1.0 to 0.0 as frequency decreases and approaches the 
critical frequency. Fortunately, periods which are free of 
scintillation do occur, particularly during years of low 
solar activity. Such periods may last for weeks or even 
months. Maximum use must be made of these favorable 


conditions when they are present. 


1.6 Terrestrial Interference 

For radio telescopes operating at below 3 or 4 times 
De interference from radio transmissions is a continuous 
problem. There are many paths which interfering signals can 
take, including single hops, multiple hops with intermediate 
reflections from the earth’s surface (possibly at almost any 
angle), or multiple hops with intermediate reflections from 
the tops of patches of sporadic E ionization. 

Often a signal will reach a receiving antenna 
simultaneously via a number of slightly different paths. The 
result is fading [24], as the signals will generally be 
different in phase. If all received components are 
approximately equal in amplitude but have random relative 
phase, the probability density of the resultant 
instantaneous amplitude A will have a Rayleigh distribution 


f(A) given by 
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Cy = SSyexp (-A2/A 2) (r1PA53) 
m ™m 


where A? is the mean square value of A. A Rayleigh 
distribution is found to give a good approximation to the 
short-term distribution of wave amplitude at long distances 
when only sky waves (ionospheric reflections) are involved 
and no ground wave is received. In the case of a strong 
undisturbed component (such as a ground wave or specular ly 
reflected sky wave) plus weak randomly scattered sky waves 
the resultant amplitude will have a Rice distribution given 


by 


ok ~(A2+B2) AB 
roa) = B exp PASEO] 1, [as (1.8) 
where ¥ = total power in signal 
A = amplitude of resultant 
B = amplitude of undisturbed component 


I. = modified Bessel function of the first kind of 
order zero 
If the random components are very small relative to the 
undisturbed component, the Rice distribution converges to a 
normal distribution. 

The above distributions are commonly observed for 
signals for which fading is present due to interference 
among multiple wavelets. The rate of such fading is 
generally quite rapid, occuring on a time scale of seconds 


to a few minutes. 
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Slower, long-term variations in signal strength occur 


due to changes in ionospheric conditions. Long-term 


variations in signal power p are experimentally found to fit 


a log-normal distribution [25] given by 


where m = 


O -& 


If signal 


= o 2 
£09) = <pe exp [Heaton fie) 


median signal power, m>0 
standard deviation of the natural logarithm of 
signal power 


strength is expressed logarithmically (in dB) then 


a normal distribution results. 


where p = 
m = 


q = 


ee Na 
rp) = Se ew PS] (ele) 


Visis 


signal strength in dB 
median signal strength in dB 


standard deviation of signal strength in dB 


Measurements of hourly median signal strengths for a given 


time of day and a given season of the year have been found 


to fit the above distribution reasonably well [26]. 


Decametric interference is difficult to escape 


entirely. 


As mentioned before, propagation via ionospheric 


reflections is generally possible only if the transmitting 


Teequencyars ‘nourmornemthane3 ori4atimestiey Most 


interference therefore subsides if the ionosphere’s critical 


frequency drops, as radio users are forced to switch to 


lower frequencies to maintain reliable ionospheric skip. 


However, ground wave communication over shorter distances is 
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still possible, so not all higher frequency transmissions 
will cease. Also, if observations are carried out at night 
when electron densities are low there is a possibility that 
the critical frequency on the other side of the earth will 
be high enough to support propagation. Signals from the 
sunlit side may find ways to reach a telescope at night. A 
prime candidate to allow such an occurence is sporadic E—E 
jonization [16]. At middle latitudes, sporadic E is more 
prevalent in summer than in winter and during the day rather 
than at night. However, it can possibly occur at any given 
time. One suspected source of sporadic E is meteor trails. 

In any case, interference to low frequency astronomy 
can potentially occur at any time and any place on the 
surface of the earth. 

Suggestions have been made that low frequency 
telescopes could be freed of ionospheric limitations by 
being built in space. An orbiting telescope, however, could 
be subject to as much or likely even more terrestrial 
interference than ground-based telescopes. One location 
which would perhaps be the most ideal is the far side of the 
moon [27]. However, such a telescope would appear to be a 
fair distance in the future. 

Another interesting idea which has been proposed is the 
creation of a temporary window in the ionosphere through 
which decametric observations could be conducted [28]. Such 
a window might be created by injecting hydrogen into the 


night-time ionosphere to reduce electron densities. 
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Interference removal during observations through an 
artificial window would be highly desirable, as interference 
could be reflected from the remainder of the ionosphere 
Surrounding the window. If one went to the trouble (and 
expense) of creating an artificial window, one would want to 


maximize the chances of successful observing. 


1.7 Previous Attempts to Deal with Interference 

Most low frequency astronomers in the past have simply 
waited for opportunities when ionospheric electron densities 
are low enough to prevent significant interference problems. 
Low frequency astronomy is most successful during a period 
of sunspot minimum. Winter nights are particularly 
attractive. However, even at such favorable times observing 
conditions are often unsuitable. For example, Dewdney [2] 
TeUPomUndueOUle Ole Cutt Ont SO sODSeCRVInNGsCUn Ing a solar 
minimum only about 12 produced results of adequate quality 
for making synthesis maps. The remaining 83% of observations 
were unsatisfactory due to interference, scintillation and 
Getrag te On. 

A few limited attempts have been made in the past to 
eliminate interference received by low frequency telescopes. 
In 1958, Shain [29] used a series of manually tuned 

Pulversei mM. cOniuUnCtion withya Na lls cross telescope .at 
1S efeavinz. Four 4.5 KHz bandwidth: filters. could be 
adjusted to interference-free gaps in a 100 KHz wide section 


Op mune spectrum cis caeNinz. Adjustment was done by an 
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operator who monitored the filter outputs on a loudspeaker. 
Shain reports that there were sentra hy a few gaps in the 
spectrum which would remain clear of interference, and the 
filter system was successful in reducing problems with 
terrestrial transmissions. 

In Tasmania ([7], [30]) a swept filter technique has 
been used to avoid interference at frequencies of 4.7 and 
10.02 MHz. A narrow filter of 2 to 3 KHz bandwidth was swept 
through a region of the spectrum 10 to 12 KHz wide at a rate 
of about 5 times per second. A minimum detection circuit 
then measured cosmic noise in clear channels between 
transmitting stations. Although this method does work, it 
limits the total observing bandwidth to a few kiloHertz and 
is relatively insensitive to interference as the spectral 
resolution is poor. Effective integration times are also 
short. 

Russian astronomers using the UTR-1 telescope at 
Grakovo [6] in 1966-68 have also reported on their methods 
of handling interference. Their receivers employed a 
variable bandwidth of from 3 to 14 KHz. The bandwidth was 
reduced at times when interference was being received. In 
addition, they mention "fast radiometer retuning’ to avoid 
terrestrial signals, evidently implying the presence of a 
human operator who monitored received signals and Kept the 
receivers away from interference. 

Human interference detectors as at Grakovo and as 


mentioned by Shain are undoubtedly very sensitive to 
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interfering signals, particularly if auditory discrimination 
is employed. However, most humans would find 12 consecutive 
hours or so (in the middle of the night) of straining to 
hear radio signals somewhat tedious. A sensitive automatic 
system which could maximize the available interference- free 
bandwidth at any given time would probably be welcome. 

With a view to characterizing the electromagnetic 
environment from 1.5 to 6.0 MHz, Wheeler [31] conducted a 
study which involved frequency-domain interference excising. 
The eventual objective of this work was the identification 
and rejection of narrowband interference from wideband 
(spread-spectrum) communications systems, which in principle 
is the same problem faced in radio astronomy. Wheeler 
recorded radio signals on magnetic tape and analyzed the 
power spectra digitally. Identifiable interfering signals 
were found to fit log-normal distributions, with 
distribution parameters depending upon the time of day and 


the season. 


1.8 Interference Detection via Spectral Analysis 

The proposed interference rejection system is based 
upon the correlation receiver, as this type of receiver is 
commonly used throughout radio astronomy for interferometry 
and aperture synthesis. 

The simplest method of detecting interfering signals is 
through spectral analysis, as a priori Knowledge of the 


frequencies and modulation of the signals is not required. 
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Terrestrial communications signals at decametric wavelengths 
are narrowband (except for military spread-spectrum systems) 
and may be distinguished from astronomical signals which are 
essentially Gaussian noise on that basis. Typical bandwidths 
of communications signals are 6 KHz for amplitude modulation 
(AM), 3 kHz for single-sideband (SSB), 1.2 kHz for 
frequency-shift-keyed telegraphy (FSK) and 100 Hz for 
continuous-wave telegraphy (CW). International agreements 
and government licensing restrict various classes of users 
to certain allocated frequency bands. 

For a correlation receiver, spectral analysis is 
possible through calculation of the cross power spectrum. 
SUCH -d spectrum consists of real and imaginary components at 
each frequency, with the real components (the in-phase or 
co spectrum) showing correlated signals which are in-phase 
Const Bi oon te and the imaginary components (the 
Quadrature spectrum) showing correlated components which 
differ in phase by 930 degrees (or 270 degrees). 

In correlation analysis the process of averaging is 
used to improve signal-to-noise ratio if uncorrelated noise 
is present and signals are relatively stationary over the 
averaging interval. In astronomy, averaging or ‘integration’ 
is usually essential as signals from sources are often very 
weak compared to uncorrelated noise received from other 
regions of the sky in, the antenna beam at low frequencies or 
compared to receiver noise at high frequencies. 


For the purpose of interference detection, the 
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sensitivity of a system employing cross spectral analysis 
will depend upon the resolution and the amount of averaging. 
For maximum sensitivity the resolution should be less than 
the bandwidth of the most distinguishing features of the 
spectra of interfering signals. AM and CW both involve the 
modulation of a continuous carrier. A steady carrier has 
zero bandwidth, implying that the sensitivity for detection 
of such signals will increase without limit as resolution is 
made finer. An FSK signal may or may not exhibit sharp peaks 
in its spectrum. A SSB signal generally will not contain any 
peaks. If peaks are not present then there is no advantage 
in having a resolution which is less than the signal’s 
bandwidth. 

Through averaging, the sensitivity will increase in 
proportion to the square root of the length of the averaging 
interval. The maximum averaging interval which may be used 
is limited by the duration of time over which the signal 
remains stationary (if the signal disappears or reverses its 
phase, continued averaging is counterproductive to its 


detection). 


1.9 A Method of Rejecting Interference 

The subject of the remainder of this thesis is an 
experimental system employing cross spectral analysis which 
was designed and built to investigate the possibility of 
detecting and removing terrestrial interference from signals 


received by a low frequency radio telescope. Spectral 
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analysis is accomplished in real time by a digital fast 
Fourier transform processor. Automatic computer analysis of 
the spectra is used to identify interference and to estimate 
the levels of broadband astronomical signals after 


interference has been excluded. 
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2. The Fast Fourier Transform Processor 


2.1 Introduction 

Spectral analysis in radio astronomy was initially 
accomplished with banks of analog filters dividing the 
spectrum into many separate channels. Analog filters, 
however, exhibit a number of problems. Precise filters 
require careful construction and are expensive. Temperature 
sensitivity and component aging cause the performance of 
analog circuits to change with time. These problems made 
astronomers look to stable digital circuits as an 
alternative. 

Digital spectral analysis was pioneered by Weinreb [32] 
in 1963. Weinreb’s technique involved the quantization of an 
analog signal into a one-bit digital signal. The digital 
signal was then correlated with delayed versions of itself 
to produce an autocorrelation function. A Fourier transform 
of the autocorrelation function produced the power spectrum. 
The advantage of Weinreb’s method is that inexpensive and 
stable digital hardware replaces expensive and unstable 
analog hardware. The method can be easily adapted to cross 
spectral analysis and its performance can be improved by 
using more levels of quantization, as for example in the 
Dominion Radio Astrophysical Observatory’s synthesis 
telescope [33]. | 

Digital correlation spectrometers of this Kind are well 


suited to high frequency radio astronomy where bandwidths of 
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many megaHertz are necessary and can be obtained by 
employing a large number of correlators operating in 
parallel. In addition, spectral lines being observed are 
relatively weak in comparison to the total power being 
received, hence the dynamic range of the correlators need 
not be large and quantization can be relatively coarse. 
Coarse quantization allows the correlators to be simple and 
inexpensive, thereby making correlation spectrometers 
economically attractive. 

In more recent years, the advent of the fast Fourier 
transform and improvements in digital technology have made 
spectral analysis by direct Fourier transformation a 
Boscom yn 4s 1S well Known, he FRivgreduces “the number -of 
arithmetic operations needed to calculate a discrete 
POUL Mer merenstorm «ot N DOInMLsS trom the order of IN@ to the 
order of N log N. FORM largea\. Savings in acompulal ton atime 
are quite considerable. 

BY reducina “the amount of computation required, the FFI 
RoSmoMenDOleN tla ys1or tae construction of simpler 
spectrometers than the correlation spectrometers described 
above. A number of FFT Spectrometers such as one at the 
Dudley Observatory [34] are already in use for radio 
astronomy. However, though the number of arithmetic 
Operations js less,’ the complexity of the operations is much 
greater. For example, the number of bits of accuracy must be 
Targer to avoid the accumulation of roundoff errors. Also, 


the control circuitry needed is far more complex than for a 
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correlation spectrometer where all stages are essentially 
identical. Thus, the greater overhead in circuitry required 
for FFT spectrometers has tended to reduce their 
attractiveness in high frequency astromomy where maximum 
bandwidth is essential. 

For interference removal at low frequencies, in 
contrast, spectrometer requirements are reversed from those 
above. As radio sources tend to be bright at low 
frequencies, wide observing bandwidth is not needed for 
adequate sensitivity. As bandwidth after a certain point 
must be obtained by adding more parallel hardware, the cost 
of spectrometers at the limits of a particular digital 
technology is essentially proportional to the bandwidth. The 
cost of a spectrometer with 50 KHz bandwidth should 
therefore be 1/100 of one with 5 MHz bandwidth. 

Wide dynamic range is essential for interference 
removal. Interfering signals can be received at levels far 
higher than background noise. Therefore, quantization of the 
incoming signals must be relatively fine to preserve 
information over the entire range of received signals. Also, 
the number of bits of accuracy retained in arithmetic 
calculations must be sufficient to not cause degradation of 
signal-to-noise ratio. 

The requirements of low bandwidth and wide dynamic 
range give the FFT spectrometer an advantage in cost and 
complexity over the correlation spectrometer for low 


frequency interference removal. 
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The fast Fourier transform emerged from obscurity in 
1965 with the publication of a paper by Cooley and Tukey 
[35]. Since that time the FFT has found widespread 
application in all areas of data and signal processing. The 
derivation and properties of the FFT are well Known and are 
Secsemibcceimemany capers samcdetexcs such as iso) =[3st. [39]. 
Ces i 142 sand [43]. Some relevant aspects of the FFI 


are discussed below. 


DeZeisealculation of the Cross and Auto Spectra 

A method will be described which allows both the cross 
and auto spectra of two real series to be calculated Wiel 
one FFT operation. Let f(n) and g(n) for OsnsN-1 be two real 
series corresponding to digitized samples of radio signals 
from two antennas forming an interferometer. Assume f(n) and 
g(n) have discrete Fourier transforms (abbreviated DFT) 


Guventpy fim) sand Gim) tor OsmsN=1. 


N-1 
F(m) = n=0 f(n) exp(-4 21anm/N) 
N-1 
E(m)x= nz g(n) exp(-j 27nm/N) eae 


Note that Fl(m) and G(m) are complex. 
The two real series f(n) and gin) may be combined into 


one complex series x(n) as follows: 


dca) = Aen) Se a SG) 
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Then by the linearity property, the DFT X(m) of x(n} will be 


X(m) = F(m) + j G(m) (amas 


Symmetry properties for the DFT’s of purely real or 
imaginary series may be used to find F(m) and G(m) from 


X(m). The symmetry property for real series states that 


F(N-m) = F*(m) 


G(N-m) = G*(m) (Zee 


Therefore 


X(N-m) = F(N-m) + j G(N-in) 


iH} 


F*(m) + 3 G*(m) 


x*(N-m) = F(m) - j G(m) (eee 


Adding equations 2.4 and 2.2 and solving for F(m) gives 


F(m) = [X(m) + X*(N-m) ] 2-85) 


No |Re 


SiiMelcmiy ssuDbtractyngsgives 


G(m) = -j/2 [X(m) - X*(N-m)] Pouce 


From the above results, expressions for the cross and 
AUMONSDeC Tra sare sreadilystound. [ifst. of all, express Atm} 


in terms of its real and imaginary components 


X(m) = X_(m) are X, (m) 


The cross power spectrum of f(n) and g(n) is given by 
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FOG) = F(m) GaGn) 


= (5 (Xm) + xem) ]1) (5/2 [xqn) - x* Orn) 1)" 
= ( = (x(m) + x*(-m) 1) (4/2 [x*(m) - X(N-m)]) 
Ey eetx(m)x (m) = X(N-myx (Nem) 


+ 4/4 [-X(m)X(N-m) + X*(m)X*(N-m) ] 


3/4 {Xr2(m) + X42(m) - X-2(N-m) - X42(N-m)} 


i] 


+ 4/4 {-[X,(m) + j Xz(m)] [X,(N-m) + 7 X; (N-m)] 


+ [X;(m) - j Xq(m)]. [X;(N-m) - 4 Xz (N-m)} 


which after some algebraic manipulation produces 


1 
a 


FG™(m) = 5 [Xp(m)Xq(N-m) + Xj (m)X, (N-m) ] 


+ 4/4 [Xp2(m) + Xq2(m) - Xy2(N-m) - X42(N-m) J CIE) 


The real components of FG*(m) are the in-phase or co 
spectrum whereas the imaginary components are the quadrature 
spectrum. 


The auto spectrum of f(n) is given by 


FF*(m) = F(m) F*(m) 
= = [X(m) + X*(N-m)] [X*(m) + X(N-m) ] 
[X(m)X*(m) + X (N-m) X* (N-m) 


X(m)X(N-m) + X*(m)X"(N-m) J 


Bile + Flr ik 


ee , 2 zs NaN RON 
= TE eee GUY Ae IG) Lar SON Ei (N-m) 


Pat Xm) oejees (in) el Xe Nom) oped <4 (Nom) | 


+ [X.(m) - 4 Xy(m)] [X,(N-m) -j Xj (N-m)]} 


which after algebraic manipulation gives 
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FF"(m) = - {X,2(m) + X;7(m) + X_° (N-m) + X; 7 (N-m) 
+ 2X,(m)X,(N-m) - 2X; (m) Xj (N-m) } (2 
Similarly, 


cS (Ge 2 .2 2(N- ;2(N- 
i xs (m) + X;“(m) + Xy“(N-m) + X4“(N-m) 


" 


- 2X,(m)X,(N-m) + 2X, (m) Xq(N-m) } (2.9) 


The above equations demonstrate the algebra used by the 
FFT processor to simultaneously calculate the cross and auto 


spectra of two real series with one FFT operation. 


2.2.2 Leakage and Windowing 

A problem which occurs during use of the DFT for 
spectral analysis is the phenomenon of leakage. An excellent 
discussion of leakage and the measures called windowing or 
smoothing used to correct it is contained in Harris [44]. 

The DFT may be viewed as a spectral decomposition in an 
N-dimensional orthogonal vector space. The basis vectors are 
of course N/2 equally spaced sinusoidal and cosinusoidal 
functions. Leakage occurs because there are a finite number 
of basis vectors whereas a physical signal, even though 
bandlimited, may have an infinite number of spectral 
components. Those components which do not precisely match 
any of the basis vectors (or equivalently are not precisely 
periodic in the finite observation interval employed) will 


exhibit non-zero projections on all of the basis vectors of 
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the DFT. The non-zero projections on DFT spectral components 
which may be far from the SReeMEREy of the actual signal are 
called leakage. 

Windowing may be used to reduce the undesirable effects 
of spectral leakage. Data prior to DFT processing is 
multiplied by a weighting function, or window function, 
which goes smoothly from zero at the end points of the data 
to 1.0 in the center of the data. The effect of windowing is 
to reduce the order of discontinuities at the boundaries of 
the data, as from one viewpoint it is discontinuities 
between the periodic extensions of the data which give rise 
to leakage. 

There are a large number of different window functions 
from which to choose, all having somewhat different 
characteristics and different effects on signals being 
processed. Harris describes many of these windows and gives 
a number of figures of merit for them which are useful for 
comparisons. One of Harris’ conclusions is that a window 
called the Kaiser-Bessel is a Superior choice for tone 
detection using the DFT. The Fourier transform of this 
window has a highly concentrated central lobe and very low 
sidelobe levels. For the problem of removing interference 
from the radio spectrum, it is crucial that a maximum amount 
of interference energy be concentrated into a few points and 
a minimum remain in the sidelobes. For this reason, the 
Kaiser-Bessel window was chosen for use with the 


inter ference-excising correlator. 
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The coefficients for the Kaiser-Bessel window are 


defined by 


w(n) = Ip{mav1.0 - (2n/N)2] , O<|n|<n/2 
To{ra] (2.10) 


where I, is the modified Bessel function given by 


to) = yg [Seed] * (2.11) 


The choice of the parameter a allows a trade off between 
sidelobe level and main-lobe width. An intermediate value of 
a=2.5 was chosen for this project. The highest window 
sidelobe level is then -57 dB and the 3 dB bandwidth of the 
main lobe is 1.57 bins, where one bin is equal to the 
difference in frequency between DFT spectral components. 

For the FFT processor, the window coefficients w(n) had 
to be quantized with 8 bits of accuracy. There was concern 
that quantization would affect the sidelobe levels. A 
Fourier transform of the quantized version of the window 
showed that there was indeed some effect, but it was not 
serious. For an unquantized window the sidelobe levels 
decrease at a rate of -6 dB per octave with frequency. The 
quantized window, however, exhibited sidelobes which did not 
decrease monotonically with frequency, but stayed relatively 
constant and varied from -70 dB to -80 dB relative to the 
main lobe level. These sidelobes are one of the factors 


limiting the dynamic range of the FFT processor. 
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2.2.3 The Effects of Windowing 

Windowing has a number of effects on signals being 
processed. Three parameters describing the most important of 
these effects are coherent gain, equivalent noise bandwidth, 
and equivalent integration time. 

The coherent gain (CG) is the gain for a purely 
sinusoidal signal, and is equal to the de gain or simply the 
sum of the window terms. Generally, CG is expressed relative 
to N, the gain of a rectangular window {w(n)=1.0 for all n}. 


Thus 
N-1 
CG = 1/N EZ) w(n) De, 


Because the width of the main lobe of the Fourier 
transform of a window function is generally larger than that 
of a rectangular window, another effect is apparent when 
broadband noise is analyzed. The amplitude of a given 
spectral estimate contains contributions from neighboring 
spectral components, resulting in an increase in the 
equivalent noise bandwidth (ENBW) of each estimate. ENBW may 
be defined as the width of a rectangular filter, with the 
same coherent gain, which would accumulate the same noise 
power as the windowed spectral estimate. An expression for 


ENBW is found to be [44] 
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For a Kaiser-Bessel window with a=2.5, the coherent 
gain is 0.44 and the equivalent noise bandwidth is 1.65. The 
major consequence of having an ENBW greater than 1.0 is a 
reduction in the signal-to-noise ratio for the detection of 
a tone in the presence of noise. For the above window the 
change in SNR relative to a rectangular window with ENBW=1.0 
is -10 log 1.65=-2.17 dB. Therefore, in order to have low 
sidelobe levels one must pay a small price with poorer SNR. 

Another result of using windowing is an effective 
reduction of integration time and a consequent decrease in 
the sensitivity of the FFI spectrometer. The decrease in 
sensitivity occurs because points near the ends of the data 
are not weighted as heavily as points near the center. Some 
of the information in the sample is therefore discarded. 

For the analysis of Gaussian noise, the result is a 
decrease in the stability (i.e. an increase in the variance) 
for estimates of noise level. The increased variance does 
not occur for individual spectral components but rather for 
the average of a number of components. Windowing causes a 
certain amount of correlation between spectral components 
which are near one another. As components are then not 
independent (as they would be without windowing) the 
variance of the average of N components is greater than 1/N 
times the variance of each individual component. 

The correlation between components of spectra of white 
Gaussian noise due to windowing has been calculated by 


Persson [45] and by Durrani [46]. Their results indicate 
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that the correlation coefficient o(k, -k, ) between power 


spectrum components separated by aon bins is given by 


pCa -k>) = w2(n) cos [2mn(kj-k2)/N| 4 


r 
| 

| Nal 

n£0 we (n) 


(2.14) 


The above result holds as long as neither component is near 
the frequencies 0 or N/2. For a Kaiser-Bessel window with 
a=2.5, the correlation coefficients evaluated as above are 
given in Table 2.1. 

An original derivation of the effect which correlation 
has upon the spectral average when used as an estimate of 
noise level is presented below. Consider an average Y of K 
components Y(k), JUsksJU+K each with variance o,? and with a 
correlation p(k, ,k, )=p(k, -k,) as above between any two 
components Y(k,) and Y(k,). Assume uU>0 and J+K<N/2 to avoid 


problems with the correlation coefficients at 0 and N/2. 


ak J+K 
Y = (1/k),,£, Y (k) 

and the variance of Y will be 
1 


soa J+K 
vi¥] = (1/K?) pBy WINGO] + BEE Covi), ¥CkQ)] 


where the double sum is over all pairs Ce with Ce ee 


The covariance is 
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Table 2.1. Spectral Correlation Coefficients for 
Kaiser-Bessel Window 


Lag Value Unset | o(k) 
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Only a few terms in the summation are required, as p(k) for 


most windows goes to zero rapidly as k increases. 


For the FFT processor in this project an average of 108 


spectral components is employed. With a rectangular window 
and white Gaussian input, all spectral components are 


uncorrelated and hence o(k)=0 for k21. 


ee ee 2 
V[Y] eee 
Rect. aR K 108 


With the Kaiser-Bessel window and p(k) as given Table 2.1, 
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The variance of Y is increased by a factor of 2.24 when the 
window is used. An equivalent integration time (EIT) may be 
defined as the reduction in integration time producing the 

same increase in the variance of the power spectrum average 
for white Gaussian noise as the use of a particular window 

function does. For a Kaiser-Bessel window with o=2.5, 


EIT=1/2.24=0.446. 
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An alternate derivation of EIT, based upon the 
time-domain correlation of two windowed Gaussian series and 
producing identical results to the derivation above, is 
presented in Appendix 1. It is found that EIT is given by a 


simple function of the window coefficients. 


a 
EIT =|n£o * (n) 
— N= 


ye) 
n=0 (ANE 
The loss of information due to windowing could be 
overcome through the use of overlapped processing as 
suggested by Welch [47]. However, overlapping increases the 
number of transforms required and therefore necessitates 
either a faster FFI processor or a reduction in sampling 
rate and bandwidth. As the decrease in variance with 
over lapping would not be as large as the increase in 
variance due to reduced bandwidth, the most efficient use of 


the set processor is with nonover lapped processing. 


2.2.4 FFT Quantization Errors 

The finite word lengths necessary in digital 
computations generally introduce quantization errors during 
signal processing. In the case of the FFT, there are two 
sources of error: coefficient quantization errors due to 
inexact representation of the sine and cosine basis vectors, 
and arithmetic roundoff errors produced by rounding after 


additions or multiplications. 
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Poe COC Tm C1relcrGuanlizalLiOn Errors 

The effect of quantization errors in the sine ais 
cosine FFI basis vectors is to convolve the Fourier 
PechicwOcieO, mtme minbUt Signal with a pattern or spurious 
sidelobes. The spurious sidelobes can be found by taking the 
Fourier transform of the quantized version of the basis 
vectors. James [48] and Tufts, et al, [49] have investigated 
spurious sidelobes in FFT’s. 

[nedacitvon «tO the mumber oT bits used “for 
coefficients, James noted that the spurious sidelobe levels 
are dependent upon whether or not coefficients identically 
equal to 1.0 are represented exactly. The use of fixed point 
binary coefficients generally results in a truncated version 
Smee eeieseprop lem IScanalogous tothe tact "that 8 bit 
binary twos complement numbers may be used for all integers 
between -128 and +127, but cannot represent +128. Improved 
performance in an FFT processor can be obtained if special 
measures are taken to recognize the coefficient 1.0 and 
ensure its exact representation. 

POMP ceo tere l cOertrcremts (/ magnitude bits plus: 1 
sign bit) the largest spurious sidelobe found by James for a 
Sorc iOU esl Coals ondbebpelow Lhermain obe., The 
levels of these sidelobes are therefore somewhat higher 
tncametnese Test! LimGemr+om a dUaltized window function as 
Gescmibec 1 SeCt 1Onweeeme. Ine side lobes-cdue torrid 
€Coermicient GUdnlizdtiom errors are a major factor in 


Minsenoecne CynaniCarendenof tnevwrrl processor. 
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2.2.4.2 Arithmetic Roundoff Errors 

A second source of error in digital FFT’ s is roundof f 
following arithmetic operations. Roundoff is necessitated by 
the finite numbers of bits available for storage in digital 
machines. Analyses of roundoff errors in references such as 
POO Feaoree S21 losis e54ised Soe anda 43iahavesshown :that 
the combined effect of many roundoff errors during an FFT 
can be represented by an addition of white noise to the FFT 
results. The level of this roundoff noise depends upon a 
number of factors including how often rounding is performed, 
the number of bits used in digital representations, and the 
accuracy of the FFT coefficients. Another important factor 
is the type of arithmetic employed. Floating point 
arithmetic allows a wider dynamic range for digital signals 
and therefore generally produces less roundoff noise than 
fixed point arithmetic where no exponents are used. An 
intermediate form of arithmetic called block floating point 
uses a common exponent for blocks of numbers and produces an 
amount of roundoff noise between that of fixed and floating 
point operations. 

With fixed point arithmetic, Oppenheim and Weinstein 
[52] show that roundoff noise is constant whereas with 
floating point or block floating point the noise level 
increases as signal levels increase. For floating point 
arithmetic and a Gaussian input signal, roundoff noise 
variance is directly proportional to input signal variance. 


Roundoff noise levels are very sensitive to the 
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rounding procedure employed. Simple truncation, for 
instance, produces high noise levels. Much better results 
are obtained with conventional rounding where if a 
fractional part is greater than 1/2 the number is rounded 
upwards or if less than 1/2, downwards. However, Weinstein 
[51] noted that the case where the fractional part is 
exactly equal to 1/2 is quite important. If this case is 
always rounded in one direction, a slight correlation is 
introduced between the roundoff error and the sign of the 
number. This correlation contradicts the assumption 
generally made in analyzing roundoff errors that the noise 
and signal are independent, and is enough to cause higher 
than expected noise levels. A "randomized" rounding 
procedure which randomly rounds this intermediate case 
upwards or downwards with equal probablilty corrects the 
situation. 

FFT roundoff noise places a lower limit on the dynamic 


range of an FFI processor. 


2.2.5 Arithmetic Overflow 
A problem which occurs during fixed point digital 
signal processing is arithmetic overflow. The addition of 
two large fixed point numbers may produce a sum which 
exceeds the maximum representable by a given number of bits. 
If it is allowed to happen, overflow causes large errors. 
Overflow is not generally a problem with floating point 


numbers because of the wide dynamic range possible. With 
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fixed point processing, however, precautions must be taken. 
The simplest precaution is to never allow numbers to grow 
large enough to produce overflow. A common method is to 
divide the results of all additions by 2. An FFT of 2” 
samples is calculated in L separate stages. In a fixed point 
FFT processor, the results of each stage of computation may 
be scaled by 2 to ensure that overflows do not occur. 

However, scaling by 2 at each stage may not be 
necessary, particularly if large sinusoidal signals which 
would cause some FFT components to grow during computation 
are not present. A better approach to fixed point processing 
is to scale only when necessary. The result is block 
floating point processing, where all computations are scaled 
an equal number of times. The number of scalings forms a 
common exponent for all FFT results. The reduction in the 
number of scalings compared with fixed point processing 
produces more accurate results and lower levels of roundoff 
noise. 

Most block floating point methods work by scaling the 
results of a particular stage only if an actual overflow is 
detected. A problem with this method is that a number of 
unscaled computations may be completed and the results 
stored in memory before the overflow occurs and the 
necessity of scaling is realized. As the previously computed 
results for the stage have been rounded before storage, a 
scaling by 2 followed by a second rounding will introduce 


additional error and lead to higher levels of roundoff 
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noise. To correct this problem, results would have to be 
stored with some of the fractional bits intact in order for 
correct rounding to be accomplished later on. In addition to 
the extra memory required for this storage, the FFT 
processor must remember which results have been scaled and 
which have not. Considerable extra hardware is therefore 
required. 

A simpler implementation of block floating point 
developed by the author involves the anticipation of 
possible overflows prior to the beginning of a stage of 
computation. For instance, if some of the results of the 
previous stage had an absolute value greater than one-half 
the maximum allowable then the addition of two such numbers 
could produce an overflow. If the possibility of overflow 
exists, all results for the stage could be scaled and the 
overflow avoided. If overflow would not have occurred then 
all new results will be less than one-half the maximum and 
scaling will not be necessary for the following stage. Such 
a block floating point scheme is less difficult to implement 
than one based on the actual occurence of overflows and 


avoids the problems of rounding described above. 


2.2.6 Overflow Correction 
The basic equation employed during FFT computation has 


the form 
Cy = Ay + ByWy - ByWy 
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as described in Section 2.3.1. The coefficients We and We 
represent points of cosine and sine enone respectively, 
and have a maximum combined value when |W [=|W, |=0.707. If 
it is assumed that the data samples AB. and B. all have 
maximum fixed point values of 1.0, then C.. can have a 


maximum value of 


C eee Orctan 0, 07) (ls ©) mm a Oe 07) CL. O)e= 20404 


Tmax 


In such a case, scaling by 2 would not be sufficient to 
prevent overflow. Even with scaling, overflow during an FFT 
is a possibility. However, it will occur very infrequently. 

The most drastic consequence of overflow is a reversal 
of the sign of the number. If overflows can be detected and 
the maximum allowable number with the correct sign 
substituted then the errors caused by overflows can be 
ecmeed appreciably. In signal processing where overflows 
may occur but with low probability, such an overflow 


correction technique performs well in controlling errors. 


2.3 Design of the FFT Processor 

Design and construction of the fast Fourier transform 
processor occupied a large part of the time spent on the 
thesis project. When the project began in 1976 there were a 
few commerical FFT processors being marketed. None of these 
machines were capable of real time processing with a 
bandwidth of more than about 20 KHz, and very few had cross 


. spectral capability. In addition, all were extremely 
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expensive and consumed considerable amounts of power. 

The availability of a number of new signal processing 
integrated circuits made the construction of a specialized 
FFT processor appear feasible, and this course was chosen. 
The design goal was a minimum of 50 kHz real time bandwidth, 
divided into 128 channels. Each channel would then be 
50 kHz/128=390 Hz in width, which is somewhat greater than 
the minimum expected bandwidth to be observed which is about 
100 Hz for marine CW transmissions. 

For a bandwidth of 50 KHz the sampling rate must be 
100 KHz, and for 128 channels the number of samples per FFT 
must be 256. For real time operation the processor would 
have to calculate a 256 point FFT in 
256x1/100 kHz=2.56 msec. A 256 point FFT can be computed in 
log,256=8 stages, each stage consisting of 128 butterfly 
operations. The time available for each butterfly is then 
2.56 msec/(8x128)=2.5 usec. A butterfly processor consisting 
of four multipliers and six adders and operating at the 
above rate forms the heart of the FFT processor designed for 
the project. An in-place, decimation-in-time FFT algorithm 
is used. A two-section input buffer allows a new set of 256 
samples to be collected while the previous set of samples is 
Fourier transformed. Eight-bit block floating point 
arithmetic was chosen for the FFT processor as a compromise 
between a simple, low power and low cost machine and one 


with wide dynamic range. 
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2.3.1 The FFT Algorithm 
The discrete Fourier transform of a series x(n), 
O<snsN-1, is defined by 
N-1 
SG Soh sayy te (218) 
where W=exp(-j27/N). This equation may be rewritten in 


matrix form as [40] 
X= [W]) & 


where X and ® are column vectors of length N and [W] is an 
NxN matrix. 

The fast Fourier transform can be viewed as a matrix 
factorization of [W]. In general if N is the product of L 
‘prime factors, then [W] can be factored into the product of 
L simpler matrices where many of the terms become zero, +1 
Ong: 

Matrix factorization of [W] works particularly well if 
N is a power of 2. In this case, [W] can be factored into 


L=log,N matrices 


[W] = [Wo JTW, IW] ALY [W, ] 


where each row of the new matrices contains only two 
non-zero terms, one of which is unity with the other being 
wk =exp(-4j27k/N) , OSkSN-1. A pecularity of this factorization 
is that in order for the matrices [W,] to be no larger than 
NxN, either the input vector x or the output vector X cannot 


remain in its natural order. Rather, the terms must be 
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permuted in a fashion called bit-reversed order. The 
‘sequence of the terms can be found by representing the 
indices i of each term, OSiSN-1, as a binary integer and 
then reversing the order of the bits. 

There are two possible factorizations of [W] which 
result in NxN matrices as above when N=2". The 
factorizations correspond to the FFT alogrithms Known as 
decimation-in-time and decimation-in-frequency. The 
decimation-in-time algorithm was chosen for this project. 

The basic arithmetic computation for a radix-2 


decimation-in-time FFT is given by 


i ae W's 


4 WRtN/2, Lg. kp (2319) 


OQ 
i} 


D=A 
where A and B are complex input values and C and D are 
complex results. Note that 


eee = exp [7 =" (ctN/2) J 


Ml 


exp( I=") exp(-j7) 


; 
-exp(-j27k/N) = ewe 


The operation in equation 2.19 is called a butterfly 
operation from its graphical representation in Figure 2.1. 


Let 
DR aay SY RET 


B= B. + 4B; 


W =o ES ee aR ey 


Then the butterfly operation can be expressed as 
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FIGURE 2.1 Complex Butterfly Computation for 
Decimation -in-Time FFT 





FIGURE 2.2 Real Butterfly Computation for 
Decimation-in-Time FFT 
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C= Cy + jCi = (Ar + BeWy - ByWy) + (CAL + BrWi + ByWr) 


H} 
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In terms of real rather than complex operations, the 
butterfly computation can be performed with four 
multiplications and six additions as shown in Figure 2.2 


56:2 


2.3.2 FFT Addressing 

An FFT is performed by executing a sequence of 
butterfly operations as described above. It is necessary, of 
course, to supply the correct data A and B plus the correct 
coefficient ws for each operation, and to store the results 
C and D in some form of memory. This control of memory 
access is provided by a scheme of memory addressing. The 
addressing scheme outlined below is an original method 
developed by the author. 

To illustrate how addressing can be accomplished, the 
example of an 8-point FFT in Figure 2.3 [37] is given. The 
dots or nodes in Figure 2.3 represent variables stored in 
memory locations. Lines entering each node from the left 
represent additive contributions to that variable, while 
arrows with a factor wk beside them represent multiplication 


by the factor. The 8-point FFT is computed in log, 8 ras: 
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FIGURE 2.3 Butterfly Diagram for an 8- Point 
Decimation - in - Time FFT 
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stages with each stage consisting of 8/2=4 butterfly 
aweratieme. Note that the results are in bit-reversed order. 
Only 8 complex memory locations are required for the above 
FFT because the results C and D of each operation can be 
stored back into the same memory locations from which the 
data A and B were read. 

A simple pattern in the addressing of the memory 
locations can be noted from Table 2.2 which shows the binary 
representations of the addresses for each butterfly in 
Figure 2.3. The A and B addresses are different in only one 
Dae 

The column of this bit starts as the most significant 
column for the first stage and moves right one column per 
stage to finally become the least significant column. 
Disregarding this particular column, the remaining bits form 
a binary counter which is incremented by one for each 
butterfly operation. 

A simple pattern is also present for the integer k 
corresponding to the coefficient ws, as shown in Table 2.3. 
The exponent k follows the pattern of a bit-reversed binary 
count, with the counting rate doubling with each stage of 
the FFT. One method of obtaining the required coefficients 
wi dumingeineseras 1setoses Llorestnescoefficients. Anza 
read-only memory in bit-reversed order, and to address the 
memory with a variable-rate binary counter where the rate is 
determined by the number of the stage of the FFT being 


computed. 
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Table 2.2. Memory Addresses for an 8-Point 
Decimation-in-Time FFT 


Data First Second Third 
Stage Stage Stage 
Address Address Address 
A G2020 OS OnmeO OR020 
B 10.80 OF a0 OR 0 aa 
A 0 0 0 0 Gast 0 
B 190 Oma Otel ae 
A Gai (0 lame On) OG 
B 110 1430 eidthoy 
A One sy Peis & Fh 16 
B ist Tis laa 11 eal 
at 
i Counts Counts 
Alternates Alternates 


Table 2.3. Powers of W for Coefficients in 8-Bit 
Decimation-in-Time FFT 


Firsitmotage Second Stage Third Stage 


0 0 0 0 0 0 
0 0 0 0 0 
00 fe Oss 
Om) ala0 1c 


2.3.3 Layout of the FFT Processor 
A block diagram of the FFT processor appears in Figure 
2.4. The processor can be broken down into six stages: 


1) analog-to-digital conversion; 2) windowing; 3) fast 
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FIGURE 2.4 Block Diagram of the FFT Processor 
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Fourier transformation; 4) power spectrum computation: 


5) accumulation; and 6) the microcomputer. 


2.3.3.1 Analog-to-Digital Conversion 

The intermediate frequency outputs of the receivers are 
bandlimited to between 416.7 kHz and 468.8 kHz (3 dB points) 
by 8-pole Chebyshev bandpass filters. These filters provide 
sharp cutoff at the bandpass edges and an asymptotic 
roll-off in the stopband of 48 dB/octave to prevent aliasing 
of out-of-band signals. The 52.1 KHz bandpass IF output 
centered on 442.75 KHz is sampled at a rate of 104.2 kHz by 
A/D converters with a resolution of 8 bits. 

Because the frequencies being sampled are relatively 
high, the times at which samples are taken must be very 
accurate to ensure that errors due to variation in the 
sampling time are small. A measure of this variation is 
given by the aperture uncertainty specification of a 
samp le-and-hold amplifier. Aperture uncertainty is defined 
to be the difference between the maximum and minimum delays 
from the time a sample-and-hold command is given to the 
point when the output ceases to follow the input. 

The aperture uncertainty time of the sample-and-hold 
amplifiers used in this project was 2 nanoseconds. The 
maximum expected error in the sampled signal due to aperture 
uncertainty, assuming a full scale signal at 450 kHz, can be 
found by noting that the maximum rate of change of the 


signal in bits/sec will be 
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= 3.62x10®8 bits/sec. 
Hence the uncertainty will be 


(2x1079 sec) (3.62x108 bits/sec) = 07/724 bits. 


The expected error is less than 1 bit. 
The A/D converters used in the project had a conversion 
time of 2.4 usec, which is much less than the 9.6 usec 


between samples. 


2.3.3.2 Windowing 

Blocks of 256 samples from each receiver are weighted 
using a Kaiser-Bessel window function to reduce leakage. The 
windowing stage is diagrammed in Figure 2.5. 

Samples f(n) and g(n) from the A/D converters are gated 
through a multiplexer onto a bidirectional 8-bit Window Data 
Bus. From this bus, the samples can be read into the Window 
Register. This register converts the sample into a serial 
bit stream and passes it one bit at a time to an 8-bit 
serial/parallel multiplier. Here the twos complement sample 
is multiplied by an unsigned 8-bit window coefficient from 
the Window Coefficient Memory. 

The resulting 16-bit serial product goes to a serial 


adder which rounds it to 8 bits. Randomized rounding is used 
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for the case where the fractional part of the product is 
exactly one-half. The rounding direction is seeermined from 
a pseudorandom binary sequence 218-1 bits long produced by 
an 18-stage shift register with feedback through 
exclusive-OR gates, as described in [57]. 

The rounded result is clocked back into the Window 
Register and is then gated via the bus to the Input Buffer 
of the next stage. A feature of the windowing stage is that 
multiplication by the coefficient 1.0 is possible even 
though 1.0 cannot be stored as one of the coefficients in 
the memory because it requires an extra bit. Because 0.0 is 
not needed as a window coefficient, a Unity Detection 
circuit which detects all 0’s at the memory’s output is used 
to inhibit the multiplication process. The original sample 
is then transferred directly to the input buffer, in effect 
being multiplied by 1.0. 

The windowed data is stored in one half of the 
two-section Input Buffer. Samples from one receiver are 
stored as the real components of the data to be Fourier 
transformed while those from the other receiver form the 


imaginary components. 


2.3.3.3 Fast Fourier Transformation 

While a new set of 256 samples is being taken, the 
previous set in the second half of the Input Buffer is fast 
Fourier transformed by the butterfly processor. A block 


diagram is shown in Figure 2.6. Address generation circuitry 
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using the pattern in Section 2.3.2 supplies data from the 
Input Buffer and FFT sine and cosine coefficients from a 
read-only memory to the butterfly processor. Butterfly 
results are calculated as in equation 2.20 and are stored 
back into the Input Buffer locations from which the data 
were read. 

During the eighth and final stage of the FFT the 
results are stored in an Output Buffer. Addresses to the 
Output Buffer are supplied in bit-reversed order, causing 
the FFT results to be-unshuffled and to appear in the Output 
Buffer in their normal order. 

The butterfly processor performs the arithmetic 
computations for each FFT and as such forms the heart of the 
FFT processor. Its operation is detailed below. 

The Butterfly Arithmetic Unit near the center of Figure 
2.6 does a butterfly operation as shown in Figure 2.2. Four 
8-bit data words A As B and BS are read into the Data 
Registers from the Input Buffer. The Data Registers convert 
the words into serial bit streams for processing by the 
serial multipliers and adders which perform the arithmetic 
computations. The devices used are Am25LS14 multipliers and 
Am25LS15 adders from Advanced Micro Devices. 

Computation results Co oe Dy and D, are passed to a 
Rounding stage where they are rounded to 8 bits prior to 
being transferred back into the Data Registers for 
conversion to parallel format. 


Rounding is an important aspect of the butterfly 
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FIGURE 2.6 Butterfly Processor 
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processor, particularly in view of the limited accuracy of 
Se bait computations. In order to keep roundoff noise to a 
minimum, no rounding is done on intermediate results such as 
products or sums during the butterfly computation. All bits 
are preserved until the last point prior to the transfer of 
results back to the Input Buffer. The level of roundoff 
noise which results is approximately 5 dB below levels 
predicted in analyses such as [50] and [51] where all 
intermediate results are assumed to be rounded individually. 
In addition to the minimum amount of rounding above, 
randomized rounding is used for cases where a remainder is 
exactly 1/2. Randomized rounding prevents roundoff errors 
from being correlated with the signs of numbers as they are 
rounded. A 20-stage shift register with feedback via 
exclusive-OR gates provides a pseudorandom binary sequence 
229-1 bits long for the determination of rounding direction. 
Special circuitry is used to ensure an exact value for 
the FFT coefficient +1.0. The Unity Detection circuit 
monitors coefficients from the FFT Coefficient Memory. The 
coefficient 80 which normally represents a twos complement 
-1.0, is also used for +1.0. By examining the addresses 
supplied to the memory, the Unity Detection circuit can 
distinguish between -1.0 and +1.0 and instructs the 
Arithmetic Unit to perform multiplications accordingly. 
Block floating point operation is implemented using the 
overflow anticipation method. The Overflow Anticipation 


circuit monitors the absolute values of numbers as they are 
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stored in the Input Buffer by an exclusive-OR of the sign 
bit and the next most significant bit. If numbers larger 
than one half the maximum absolute value are detected, 
scaling is performed in the next stage of FFT computation. 
The Exponent Counter keeps track of the number of scalings 
so that results can be normalized later on. 

If an overflow should occur, the Overflow Detection and 
Overflow Correction circuits respond by substituting the 
maximum allowable number of the correct sign for the 


overflow result. 


2.3.3.4 Power Spectrum Computation 

The power spectrum computation stage is outlined in 
Figure 2.7. It consists of the 4-part Output Buffer, two 
8-bit busses leading to two registers Y1 and Y2 and 
connected by transmission gates, an 8-bit parallel 
multiplier, a 23-bit adder/subtractor and accumulation 
register, a storage register, and assorted control and 
addressing circuits. 

By gating the correct FFT results from the Output 
Buffer into registers Y1 and Y2 and adding or subtracting 
the resultant multiplier products appropriately in the 
accumulator, all] the required auto and cross spectrum 
components can be computed. The computations are done 
decorenndm: Omedudtaonsmcm (smc sondndecn J) excepts thats sca ling 
by the factors 1/4 and 1/2 in these equations is not done 


and the auto spectra are computed in complementary form. 
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FIGURE 2.7 Power Spectrum Computation 
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Correction for these discrepancies is made later by the 
microcomputer. 

Control signals for the power spectrum stage are 
generated with a combination of random logic circuits and 
fast read-only memories. The multiplier and accumulator 
sections are implemented using a TDC-1008u 


multiplier/accumulator from TRW Semiconductors. 


2.3.3.5 Accumulation 

The process of accumulation is equivalent to the 
integration or averaging of a large number of spectra. 
Averaging is commonly used in cross spectral analysis to 
improve the signal-to-noise ratio when weak signals are 
combined with high levels of uncorrelated noise. As each FFT 
requires 2.458 msec, the FFT processor is capable of 
averaging (60 sec/min) x (1/2.458x10-3 sec) = 24,400 spectra 
per minute. 

A block diagram of the accumulation stage is shown in 
Figure 2.8. Power spectrum results from the previous stage 
are first of all scaled in order to align decimal points 
with those of previously accumulated results. Scaling is 
necessary because the output of the FFI processor is in 
block floating point form with a common exponent. The value 
of the exponent from the FFT processor is used in the 
scaling step to normalize the power spectra. Scaling is 
accomplished by simply shifting the power spectra in a 


48-bit shift register a number of places according to the 
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exponent value. 

The speed of the scaling operation is enhanced by using 
a bidirectional shift register and loading the power spectra 
in the center of this register so that the total number of 
shifts required is never more than half the total scaling 
range. In addition the shift registers are arranged to shift 
results by two places on each clock cycle because the 
required number of shifts is always even (i.e. a shift of 
one bit position in the FFT results produces a shift of 2 
places in the power spectrum, as the power is a product). 
These measures to provide fast scaling are necessary for the 
accumulation process to Keep pace with FFT production. 

A scaled power spectrum is passed to a 24-bit adder 
which adds the new spectrum to the accumulating total of 
previous spectra. The actual length of the numbers being 
added is 48 bits which the adder handles in two steps, 
saving the carry between add cycles. The large number of 
bits used during accumulation is necessary because of the 
long accumulation times. 

Accumulation results are stored in a charge-coupled 
device (CCD) shift register memory with a capacity of 3072 
bytes. The memory can hold 128 spectral components, each of 
48 bits, for the two auto and two cross spectra. The CCD 
memory has advantages over random access memory of low power 
consumption, small size and fast read/modify/write 
capability. 

The 24-bit Register in Figure 2.8 serves to latch old 
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CCD contents for input to the adder during accumulation and 
Gioome_Omie lds accumilalion resulis during Gata transrers to 
iLnemimtenocomputer. Data transfers occur after-a specitric 
number of spectra have been added. The CCD contents are then 
zeroed and a new accumulation begins. 

Results are transferred directly to the microcomputer’ s 
memory using a direct memory access (DMA) technique. The FFT 
processor gains control of the microcomputer data and 
address busses, and supplies the addresses to write 
accumulated spectra directly into memory. The complete DMA 
transfer requires one FFT cycle, or 2.458 msec. 

The accumulation stage performs one more function, 
name ly aynetnoomeues demodulation. In order to eliminate 
crosstalk between receivers in radio astronomy, a technique 
called switching or synchronous modulation and demodulation 
is commonly employed. The FFT processor accomplishes 
switching by inverting the phase of the RF signal from one 
of the antennas with each consecutive FFT. The switching is 
equivalent to modulation by a bipolar square wave with a 
frequency of about 203 Hz. 

Modulation is accomplished using a diode ring balanced 
mixer driven by current sources of equal amplitude but 
SopOoeeesiom. (Ne umiserdnde gail Uimesyot tne ssquare wave 
pee cacer ol iy Controlled to neduce problems Of Panging in 
the receiver filters and other undesired modulation of the 
facomimeesiOndis. ohentami ng of the square wave is adjusted 


so that receiver and cable delays are compensated for, and 
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the phase transitions at the A/D converters coincide with 
the beginnings of new blocks of FFT samples. Further 
discussion of the modulation circuit is contained in Section 
3 H4., 

In order to undo the above modulation, the accumulation 
stage performs synchronous demodulation. Rather than adding 
consecutive spectra, the add is changed to a subtract on 
alternate spectra for the in-phase and quadrature spectra’. 
Crosstalk occuring after the modulation thereby tends to 


cancel itself. 


2.3.3.6 The Microcomputer 

The final stage of the FFT processor is the 
microcomputer. A block diagram of the computer system 
appears in Figure 2.9. 

The microprocessor which forms the central processor of 
the microcomputer is a Motorola MC6800. The microcomputer 
can be operated with either of two operating systems. One of 
these is Motorola’s MIKBUG which uses a teletype, and the 
other is an operating system written by the author for a 
small console with a 30-key Keyboard and 7-segment readouts. 

Programs for the microcomputer were written in assembly 
language on the University’s central computing facility. 
Machine code for the programs was then loaded into the 
microcomputer via telephone using a modem and acoustic 


coupler. An audio tape interface allowed local recording and 


‘Synchronous modulation is possible for the cross spectra 
only, as the auto spectra are always positive. 
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FIGURE 2.9 The Microcomputer 
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playback of programs and data. 

One of the features of the microcomputer is a high 
speed arithmetic processor. This processor, an Am9511 from 
Advanced Micro Devices, allowed fast 32-bit floating point 
arithmetic calculations for robust estimation after spectra 
were received from the FFT processor. The improvement in 
speed over software arithmetic was a factor of about 100. 

The microcomputer has 64 lines available through 
Peripheral Interface Adapters for parallel digital input or 
output. Half of these were required for interfacing to the 
FFT processor during normal operation. The remaining lines 
could be connected to various points in the FFT processor 
via ribbon cables during testing and debugging. 

Eight channels of digital-to-analog output were 
included. These were used to drive the X and Y axes of an 
oscilloscope to display spectra as they were received, and 


also to drive chart recorders for a record of fringes and 


other information such as the number of points being deleted 


and the amount of clipping encountered. 

A serial interface to the programmable frequency 
synthesizer used for the second local oscillator was 
provided to allow control over the receiver center 


frequency. 
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2.4 General Construction 

The FFT processor was designed with two major goals in 
mind: cost-effectiveness and low power consumption. Power 
consumption was a consideration in the early stages of the 
project because initially the FFT processor was to be 
operated at a remote observing site with no available 
commercial power. 

The processor used a combination of standard TTS 
Schottky TTL, low-power Schottky TTL, MOS and CMOS digital 
integrated circuits. Generally the type of logic consuming 
the least power was used wherever possible. Schottky TTL was 
employed in sections with fast clock speeds (up to 20 MHz) 
and where short propagation delay times were critical. 

The circuits were wire-wrapped on a number of 
approximately 10 inch by 10 inch perforated boards. 
Wire-wrapping allowed dense packing of the IC’s and also 
easy wiring changes. The circuitry was designed and tested 
in small sections. High speed sections in particular were 
kept compact to minimize the lengths of interconnections. 

Good grounding is an important consideration in digital 
design. The best approach is a solid ground plane for low 
inductance, but this would have been difficult to implement. 
A satisfactory approximation to a ground plane was the use 
of two grids consisting of interconnected heavy bus wires to 
which ground and V_. connections were made. Bypass 
Capacitors at regular intervals on the grid Kept noise to a 


minimum. 
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A problem commonly encountered in digital systems is RF 


radiation and consequent electromagnetic interference with 
communications services. In a radio astronomy environment, 
electromagnetic interference can be disastrous. Spectrum 
analyzer measurements showed that the FFT processor and 
microcomputer produced considerable RF radiation, with 
harmonics up to many hundreds of megaHertz. The highest 
levels of RF happened to be right at 22 MHz, creating the 
spectre of an interference-excising system with a paramount 
purpose of removing self-inflicted interference. 

The solution for this disturbing problem was the 
placement of the FFT processor and the microcomputer inside 
a shielded enclosure built from copper screen. Al] 
connections passing through the enclosure were low-pass 
filtered by LC sections, feedthrough capacitors and ferrite 
beads to prevent RF from escaping. The shielding and 
filtering were quite successful in eliminating 


electromagnetic interference problems. 


2.5 Laboratory Tests of the FFT Processor 

A number of laboratory tests of the FFT processor were 
conducted in order to verify the correct operation of the 
hardware and to measure the performance of the completed 


machine. 
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2.5.1 Comparisons to Computer Simulations 

A software model of the FFT processor was generated in 
FORTRAN for testing purposes. The model could simulate 
combinations of sinusoidal and Gaussian noise signals and 
calculate cross and auto spectra using arithmetic identical 
to that of the FFT hardware. The same signals were then fed 
to the FFT processor and its results compared to the model's 
results. By this method, faults in the hardware could be 
found very easily. 

The microcomputer was used as the control element for 
the FFT hardware during these comparisons. Test points were 
provided throughout the processor to allow the microcomputer 
to control hardware operation, to provide simulated data, 
and to monitor results. Major sections of the hardware could 
be tested independently or in tandem, allowing faults to be 
isolated very quickly. 

Due to the complexity of the entire machine, the above 
testing method proved to be invaluable both during the 
original construction and debugging of the hardware and for 
maintenance following construction. The software model also 
allowed aspects of the FFT processor to be tested before a 


commitment was made to hardware. 


2.5.2 Windowing 
Four photographs showing the spectrum of a sinusoida!l 
signal as analyzed by the FFT processor appear in Figures 


Dea Oe andee. 11. lm photowAyof 2. 10ma rectangular window was 
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Figure 2.10. FFT Processor Spectrum of Sinusoid Using 


Rectangular Window 
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Figure 2.11. FFT Processor Spectrum of Sinsoid Using 


Kaiser-Bessel Window 
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employed with the frequency of the sinusoid being exactly 
equal to the frequency of one of the FFT basis vectors. No 
significant sidelobes were present. In photo B, the input 
frequency was exactly between two of the FFT basis vectors 
and sidelobes appeared with the characteristic (sin x/x)2 
amp 1itudes. 

In Figure 2.11 a Kaiser-Bessel window was employed for 
the same sinusoidal signals as above. Sidelobe levels due to 
leakage were now controlled. Spurious sidelobes due to FFT 
coefficient quantization may be seen in 2.10 A and 2.11 A 
and B. The largest of these occurs at 47 KHz and is between 


56 dB and 60 dB below the peak of the signal at 5 kHz. 


2.5.3 Phase Relationships of the Co and Quadrature Spectra 

Two tests, one with a variable phase Sinusoidal source 
and one with delayed Gaussian noise, were conducted to 
verify the correct phase relationships of cross spectral 
components. 

The results of the variable phase sinusoid test are 
illustrated in Figure 2.12. The powers of the peaks of co 
and quadrature spectra are plotted against the phase 
difference between the two sinusoids correlated by the 
processor. As expected, the co and quadrature components are 
always 90 degrees apart. 

For the second test, noise from an audio noise 
generator was fed in parallel into two 512-stage 


bucket-brigade analog shift registers. The clock rates of 
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the shift registers could be controlled independently to 
produce any desired relative delay between the register 
outputs. A delay +t should produce a phase shift o which 


varies in proportion to frequency f, 
> = 2nft 2m) 


Figures 2.13 and 2.14 show the spectra produced with 

Bold GivemOelay Seon JUMUSeC sande A00RNSec, arespectivelly.. ine 
in-phase and quadrature spectra have cosinusoidal and 
sinusoidal variations of amplitude vs frequency', as 
expected from the phase difference vs frequency in equation 


agents 


2.5.4 Roundoff Noise 

Roundoff noise is one of the factors which places a 
lower limit on the dynamic range of the FFT processor. In a 
block floating point design, the level of roundoff noise is 
dependent upon the amount of scaling which occurs during 
processing and is therefore dependent upon signal levels 
(particularly sinusoidal signals as their power is 
concentrated in a few spectral components). Measurements of 
roundoff noise were made both with and without the presence 
of sinusoidal signals. 

Figure 2.15 shows a plot of equivalent input noise 
levels for the auto and cross spectra vs actual correlated 


noise input levels (no sinusoids present). As the input 


‘The scale in the figures is logarithmic, hence the peaks 
appear to be flattened. 
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FIGURE 2.13 Cross Spectra of Delayed Noise ( Delay = 100 jrsec ) 
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FIGURE 2.14 Cross Spectra of Delayed Noise (Delay = 200 sec) 
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level decreases, the auto spectra depart from a linear curve 
(dashed line) much earlier than the cross spectra. The 
reason is that roundoff noise for the auto spectra is fully 
correlated and its power has a positive mean value which 
appears as an additive contribution to the auto spectra. 

For the cross spectra, however, roundoff noise is 
uncorrelated and has zero mean. The cross spectrum therefore 
more accurately reflects the true input noise level, 
departing from it only as levels become very small. The 
curve for the cross spectra approaches a lower limit equal 
to the quantization noise of the A/D converters. As the 
maximum allowable input levels are +0.5 volts and the number 


of bits is 8, the quantization step size S is 


Seely Ol hie. Oo Lamy. 
wee 


The well Known theoretical rms value of quantization noise 
is then 
aS 
% Se ee LON cms 
¥12 
As the signals being correlated in Figure 2.15 are 
identical, the quantization noise due to A/D conversion is 
correlated and appears as an additive contribution to the 
in-phase cross spectrum. 
A number of measurements of roundoff noise in the 
presence of sinusoidal signals are plotted in Figure 2.16. 


The noise level is approximately proportional to the 


sinusoidal level. The method employed in making these 
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measurements is outlined below: 


ig. 


The signals analyzed consisted of a common sinusoid to 
which uncorrelated Gaussian noise was added in two op 
amp summers. 

After analysis by the FFT processor, the frequency 
components containing the sinusoid were deleted from the 
auto spectra and averages of the remaining components 
taken. 

The averages of the auto spectra when analyzing the same 
level of uncorrelated noise alone were subtracted from 
the averages in step 2 to find the excess noise due to 
roundoff. 

An equivalent input noise level was found which would 
produce the same changes in the auto spectra as the 
roundoff noise in step 3. 

As roundoff noise is not windowed, it exhibits no 
correlation from one component to another. Its effect on 
the variance of cross spectral estimates is therefore 
not as large as that of the windowed input noise found 
in step 4. This noise level was divided by 4/EIT to find 
an equivalent windowed input noise level which would 
cause the same increase in variance of cross spectral. 


estimates as the roundoff noise. 


The addition of extra noise is necessary because without it, 


many of the spectral components containing roundoff are so 


small that they are rounded to zero. A true measurement of 


the roundoff noise is then impossible. 








jone)pet?t ell . ws soon THA ens xd 


; . oe , pale 


“i : 7T £ og 77 a0 = ~~ 420 surwe ont gatntstnad 4 Jae HF 


vs ee 
2inerognes oetriheder oll Je “sapatove f ae B99 


vY ea 















= 
P 
‘es 





ome ant onto! sme. re end eqe olw: artt/#e 
Stew enots seton t walle 


oq aw) 32 Ton eae JAD ay T brvi T 3 iN dere nt: 
= iz 


“itriees St .Sqwotirtitye ten ‘ef set 
an | 
ne foe?%> all .setiene of taaramieo ona. 


asctsiati ef aie! tae Tentyega seni % oe 
‘ . “oom 
Eruct s2fonm -wonl Cowie ona ts tad} 4 
_ ae 
‘pett of T. 7)" yo bebiwth gem Havel saben a 


= i Laguw tat rw hawt - 2ion tence bin oF a 
: ae 


Jeviasqe 26070 to cons} a6 int sexes nt , | 
i wow are 7 aN 7 _ Pad 2e 
ape nm 


4 , | " i 
til As TC an & Ti. > | 
4 Z ay ed 


at sucrit tw eeusoed un 
ee a : a 

of 7% Vice yor pean allt 

ay we > ses Aen ; 





7 a 










: ' ov i 


A chee 


2.5.5 Correlation Tests 

A number of long term correlation tests of two 
independent noise sources were conducted to investigate the 
FFT processor’s performance as a correlator. A small offset 
was found to occur when synchronous demodulation or 
switching was not employed, but otherwise the correlations 
were as expected from the theory described in other sections 
of this thesis. . 

Correlation tests were carried out using two General 
Radio 1390B random noise generators as noise sources. The 
number of spectra accumulated per integration was K=105, 
equivalent to an integration time of 
105 x 2.46X10-% sec = 246 sec. An average of 30 trials with 
K=105 was taken making the total integration time 123 
minutes for each test. 

Equivalent levels of correlated input noise for a 
number of levels of uncorrelated input with and without 
switching are shown in Table 2.4. 

With switching no detectable correlation is present in 
the cross spectra. Without switching a small offset occurs, 
probably due to either stray pickup of a common signal by 
the input lines or a small amount of crosstalk between the 
A/D converters. The first is more likely because the offset 
is independent of level. 

The relationship between the levels of the auto spectra 
and the variance of the cross spectra for uncorrelated noise 


inputs is shown in the results in Table 2.5. As derived in 
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Section 4.2.3 the product of the means of the auto spectra 
divided by 2K gives an accurate estimate of the variance of 
the cross spectra. 

Finally, in Table 2.6 the variance of the means of 30 
cross spectra is compared with the expected variance if all 
components are independent. The average ratio of the actual 
variance to the expected variance if all components are 
independent is 2.25 which is very close to the predicted 
increase in variance due to windowing found in Section 
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Table 2.4. Correlation Test Results 


Uncorrelated Switching Equivalent 
Input Level Correlated 
(mVrms ) In-Phase 


Level (mVrms) 


100.0 YES 2841.54 
31726 ies 262.49 
10.0 TBs WOW SE 5 
100.0 NO 1.14#1.54 
Sy ol 8 NO 802.49 
10.0 NO » [fae Ae 


Equivalent 
Correlated 
Quadrature 
Level (mVrms) 
165s1e, 54 

mat 49 


044.15 


oo tn.5 4 
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Table 2.5. Auto Spectrum Levels and Cross Spectrum Variance 


(K=405) 


Input Autol Auto2 AUTO nex In-Phase 


Level Mean Mean Auto2 Variance 
(mVrms ) OK 
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Table 2.6. Variance of the Means of the Cross Spectra 
(K=105, No. of Components = 120) 


Input ] Variance of Variance of Ratio 

Level or Cross Spectrum Means of 30 

(mVrms) Q +120 Cross Spectra 
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3. Receivers 
A Dlock ‘diagram of the receiving system used for the project 
appears in Figure 3.1. The antennas, RF amplifiers, first 
mixers, first IF amplifiers and first local oscillator were 
part of a previous project at the Dominion Radio 
Astrophysical Observatory. This project was the first low 
frequency earth rotation aperture synthesis mapping of the 
north polar sky and is described in Dewdney [2]. 

The present system used the 300 KHz bandwidth outputs 
Ctimcneeriurst. dteampbivfiers, at 5.0 MHz-to feed second mere 
and second IF amplifiers, forming a sharply bandlimited 
signal from 416 KHz to 468 KHz suitable for sampling and A/D 
conversion. The receiving equipment designed and built 
during the course of the present project is described in 


Lise chdaplen. 


3.1 Second Mixers and IF Amplifiers 

A schematic of a second mixer and IF amplifier is shown 
ine GUPC. o.- 2). 

The mixer consists of an MC1496 balanced modulator. 
This device was chosen because of its good linearity and 
michevsOlation |>608dB) between) RF, IF and#iocal oscillator 
ports. The high isolation reduces crosstalk between the two 
receivers via the common local oscillator. Isolation is also 
improved through the use of a balanced hybrid splitter for 
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FIGURE 3.1 Receivers 
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The 440 KHz IF amplifier following the mixer consists 
of three LM371 cascode amplifier stages and an output buffer 
stage. Coupling between stages is accomplished via 
double-tuned tank circuits. The 8 tank circuits form an 
8-pole Chebyshev bandpass filter with a design ripple of 
O01 dBe 

With the Q of about 150 available in commercial IF 
transformers at 455 KHz, a Chebyshev filter with 0.01 dB 
ripple was the best filter characteristic attainable. Sharp 
cutoff at band edges can be obtained in filters by the 
addition of zeroes in the transfer function at the edges, as 
for example in elliptic filters [58]. However, such zeroes 
cause ripples in the out-of-band response. It was deemed 
more desirable to minimize the response to out-of-band 
signals as much as possible rather than having an extremely 
sharp cutoff. For this reason, a minimum-phase Chebyshev 
filter seemed most suitable. 

In order to achieve a symmetrical bandpass response it 
was necessary to use inductive coupling for the first two 
double-tuned circuits and capacitive coupling for the second 
two. Such a combination is necessary to get the correct 
number of zeroes in the bandpass transfer function. Very 
briefly, consider an 8-pole low-pass Chebyshev transfer 
function [59] 
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(s2+w2)/(s BW) and produces a bandpass transfer function 


with 8 zeroes. 
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The correct number of poles can be obtained by 
cascading the transfer functions of 4 double-tuned stages, 
each contributing 4 poles. However, with either transformer 
or inductive coupling a double-tuned stage contributes only 
one zero. The transfer function is 


sd SK 
G, (s) = B 


st + cuss ~ C48° Dia ome 
A cascade of four such stages would result in an overall 
G,p (s) with only 4 zeroes. The result is a somewhat skewed 
frequency response which is particularly noticeable when the 
filter bandwidth is an appreciable fraction of the center 
frequency. 
A capacitively coupled double-tuned circuit, on the 


other hand, exhibits 3 zeroes: 


5 °K 
G(s) = C 
cs) 


i 


Seed aoe AC coe ad Sak id 


By using two capacitively coupled and two inductively 
coupled stages, 8 zeroes can be obtained. 

The desired filter response was obtained by the choice 
of component vaiues in the double-tuned stages. The 
calculations are lengthy and will not be included herein. A 
short description of the design and tuning procedure 


follows. 
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Each tank circuit inductor could be tuned individually 
to allow pole locations to be adjusted. However, since no 
identical poles were needed, overcoupling of the 
double-tuned circuits was used to ensure maximum energy 
transfer. If the two tank circuits of an overcoupled 
double-tuned stage are tuned identically, then the actual 
double-tuned poles occur at different frequencies, with the 
difference being dependent upon the degree of overcoupling. 
Element values in the second IF amplifiers were chosen to 
produce differences in poles of slightly less than the 
desired values so that tuning could be used to set the poles 
exactly. 

Pole damping to obtain a Chebyshev response was 
accomplished by resistors marked with asterisks in Figure 
oleae Approximate resistor values were found through 
calculation but the final selection was made experimental ly 
by physically measuring the Q’s of each tank circuit with 
all components except the damping resistors connected. 

Interaction between the two halves of the double- tuned 
circuits did not present any problems during tuning of the 
IF amplifiers. It was possible through careful adjustment of 
the inductor cores and by some trial-and-error fine 
adjustment of resistor values to obtain a frequency response 
with less than 0.1 dB of ripple. Phase matching of the two 
IF amplifiers was to within 5 degrees. 

One problem which was encountered was a lack of 


mechanical stability in the IF transformer cores. The 
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amplifier frequency response was found to change, 
particularly in response to physical shocks or movement. 
Regular monitoring of the frequency response was necessary 
to maintain proper tuning. 

The final stage of the second IF amplifiers is an 
output buffer consisting of an LM318 operational amplifier 
and an LHOOQ2 current driver. This stage is capable of 
driving a 50 ohm load at high levels (+10 volts) with low 
distortion. The required input level for the A/D converters 


Ts smuchmless than this i(+0.5' vo tts); 


3.2 RF Amplifier, First Mixer and First IF Amplifier 

A schematic of the RF, first mixer and first IF stages 
is shown in Figure 3.3. The circuit is a modified version of 
the receivers originally used for the north polar synthesis i 
telescope. 

The original receivers used a bipolar transistor in the 
first stage of the cascode RF amplifier and a bipolar 
integrated circuit (an MC1550) as the mixer. These bipolar 
stages were found to produce high levels of intermodulation 
distortion when observing was attempted during periods of 
strong interference. 

Bipolar transistors are highly susceptible to 
intermodulation problems because their large-signal i-v 


characteristic is exponential in nature [60]: 


ve = ay, lexp(V_/V,) 2 iki) ae Azo lexp (V/V) — 1.) 
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where I, = collector current 

@511 Ao. = Constants depending upon transistor 

construction 
Vz(V~) = voltage drop across emitter (collector) 
junction 
Vr = voltage equivalent of temperature 
= temperature/11600 K/volt 

A series expansion of this characteristic contains many 
odd-order terms which give rise to intermodulation. 

Field effect transistors are well Known to produce low 
levels of intermodulation distortion. The i-v characteristic 


for a FET is a simple square law [61]: 
oe hy gee alee eal 


where Ins = saturation drain current 
Ves = gate-source voltage 


Inss = ts ly 20 

ve = pinch-off voltage 
As no terms of higher order than two are present (at least 
ideally) FET’s can amplify or mix high-level signals with 
little intermodulation. 

Replacement of the RF amplifier and mixer stages with 
FET’ s as shown in Figure 3.3 reduced problems of 
intermodulation in the receivers considerably. The third 
order intermodulation intercept point (as measured in a 


standard two-tone intermodulation test) was found to 


increase from -60 dBm to -11 dBm. 
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The noise temperatures of the two receivers after 
modification were measured to be 900°K and 1500°K. 

The first IF amplifiers consist of two cascode 
amplifiers (MC1550’s) and an output buffer stage. The 
double-tuned circuits which provide interstage coupling are 
adjusted for a Butterworth filter response with a bandwidth 
of 300 KHz centered on 5 MHz. Following front end 
modifications the responses of the two RF and first IF 
stages together were adjusted to be matched within +0.15 dB 
in amplitude and +1.25 degrees in phase over the central 
200 KHz of the pass band, as it was only this region of the 


First IF output which was actually used. 


3.3 Input Bandpass Filters 

Most of the RF selectivity of the receivers was 
provided by the tuned circuits associated with the RF 
amplifier stage in Figure 3.3. These tuned circuits were 
adjusted for a 3 dB bandwidth of 600 KHz. 

Additional out-of-band rejection was provided by 
external bandpass filters. These filters had a bandwidth of 


5 MHz centered on 22.25 MHz and were specially designed for 


high rejection at frequencies far below and above the center 


of the pass band. The filters had a 3-pole Chebyshev 


response with 0.01 dB ripple and are shown in Figure 3.4. 


The filter’s rejection at 1 MHz was greater than 100 dB 


and from 50 MHz to 350 MHz was greater than 50 dB. To obtain 


good performance over such a wide frequency range required 
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selection of components with low parasitic impedances. For 
example, the 680 pf capacitors were chip capacitors with 
extremely low inductance and the 0.074 uH coils were toroids 
with high Q and a high self-resonant frequency. High Q 
inductors limited the insertion loss at 22.25 MHz to 1.0 dB. 
Construction was also crucial. Low frequency rejection 
was found to be highly dependent on good grounding and the 
use of a solid ground plane. Good high frequency rejection 
required shielding between the two halves of the filter to 


prevent radiative coupling. 


3.4 Synchronous Modulation 

Synchronous modulation of the RF signal from the east 
antenna was acccomplished with a Hewlett Packard 10534A 
balanced modulator. The diode ring in this modulator was 
driven through a length of coaxial cable from the 
observatory building by a reversible current source 
diagrammed in Figure 3.5. 

The phase switch control input was a binary signal from 
the FFT processor which alternated between +12 V and -12 V 
ateascate sOfec0o Hz ane, .mputacircult: consisting of <an 
adjustable RC network and back-to-back Zener diodes allowed 
a variable delay to be introduced. The delay was adjusted to 
produce phase transitions of 180 degrees at the east 
receiver A/D input which coincided with the beginnings of 


new blocks of FFT samples. 
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The first op amp formed a Schmitt trigger with its 
output limited to +8.8 volts by feedback through Zener 
diodes. The second op amp was a voltage-to-current converter 
which responded to its +8.8 volt input with an output of 
+10 ma. 

An RC network at the output limited the rise and fall 
times of the current output to reduce ringing, as the output 
drives a long length of cable. 

A low-pass filter just before the balanced modulator 
prevented RF signals being conducted by the synchronous 
modulation cable and introducing crosstalk between the 


receivers. 


3.5 Local Oscillators 

The first L.G. was a crystal controlled oscillator at a 
frequency of 27.25 MHz. This oscillator signal was 
distributed from a central point via separate buffer stages 
and cables to the east and west receivers. The RF, first 
mixer and first IF stages of each receiver were physical ly 
located near their respective antennas and were connected 
via coaxial cable to the second mixers and IF amplifiers 
which were located in the observatory building. 

The second L.O. consisted of a digitally controlled 
Fluke 6039A frequency synthesizer. A serial digital 
connection to the microcomputer allowed the frequency of the 
second L.0. to be changed under computer control. As the 


synthesizer required about 40 parallel control lines, a 
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serial-to-parallel interface employing shift registers and a 
simple binary pulse width modulation scheme was built to 


reduce the number of control lines needed to only one. 
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4. Robust Estimation 


4.1 Introduction 

In this chapter, an estimation procedure well suited to 
the problem of minimizing the effects of terrestrial 
interference in the cross power spectrum is derived. Section 
4.2 begins by defining the estimation problem and finding 
the expected probability density functions of the cross and 
auto spectra. Historical aspects of rejection rules and a 
branch of statistics Known as robust estimation are 
discussed in Section 4.3, and the concept of robust 
maximum-]likelihood estimation (M-estimation) is introduced. 
The characteristics and performance of some M-estimators are 
presented in Section 4.4. A short discussion of adaptive 
estimation follows in Section 4.5. Finally, Section 4.6 
proposes a robust estimation procedure combining rejection 
rules and M-estimates which the author believes is a good 


choice for excising terrestrial interference. 


4.2 Definition of the Estimation Problem 

The underlying assumption which makes the removal of 
terrestrial interference from the cross power spectrum 
possible is that the interference will be narrow band and 
will contaminate only a certain percentage of the components 
of the spectrum. At least some of the components at any 
given time will be free of interference and these 


components, if they can be identified, can be used for 
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astronomy as if there were no interference at all. 
Components with no interference can be expected to 
exhibit a normal probability distribution, the parameters of 
which are found below for both the cross and auto spectra. 
Components containing interference will deviate from the 
normal distribution, introducing contaminated points and 
heavier tails in the distribution. Such contamination would 
be disastrous if a simple mean of all points was used as an 
estimate of the center of the distribution, as the mean is 
very sensitive to outlying points. The robust estimation 
procedures described later in this chapter are designed to 
provide protection from outlying points when estimating the 


center of a contaminated normal distribution. 


4.2.1 Probablility Density Function of Averaged Cross 
Spectra : 

This section will show that, in the absence of 
interference, an averaged cross power spectrum will have a 
normal probability density function (pdf). Consider the two 
signals being correlated as two real series f(n,k) and 
g(n,k), where OSnSN-1 and 1<k<K. The index n represents 
points within a block of FFI samples, whereas the index k is 
incremented with successive FFT’s. The FFT’s of f(n,k) and 
g(n,kK) for a particular value of kK will be called 
instantaneous FFT’s, denoted by F(m,k) and G(m,k), where 
OsmsN-1 and F and G are complex. 


The instantaneous cross power spectrum, denoted 
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* 


FG (m,kK), will be 


FG*(m,k) = F(m,k) G*(m,k) 
= Fr(m,k) Gy(m,k) + Fy(m,k) Gy (m,k) 


+ j[Fy(m,k) Gr(m,k) - Fr(m,k) Gi(m,k)] Cart} 


The subscripts r and i denote real and imaginary components. 
Before determining the pdf’s of the real and imaginary 
components above, a number of assumptions will be stated. 
First, f(n,k) and g(n,k) consist of samples of Gaussian 
noise bandlimited to exactly the Nyquist frequency. Second, 
the correlation between f(n,k) and g(n,k) is identical in 
amplitude and phase at all frequencies. Then F_(m,k) and 
F. (m,K) will be normally distributed with mean 0 and 


variance 9,2, denoted as N(0, 5,2), and G_ (m,k) and G, (m,k) 


will be N(0, o,2). The relationships between the four 


ie 
components above can be seen in the variance-covariance 


and PQ refer to the 


correlation coefficients of the real and imaginary 


matrix in Table 4.1. Here, Po 
components of the cross power spectrum (the co and 
quadrature spectra). 

From Table 4.1 it is possible to calculate the expected 
mean and variance for the cross power spectrum. Consider the 
real part first, pa th Gy Using results derived in 
Appendix 2 for the product of two Gaussian variables, the 


expected mean and variance will be 
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Table 4.1. Variance-Covariance Matrix for the Complex 
Spectra of Two Partially Correlated Gaussian 


Signals 
F F G G 
1g 1 je 1 
F 
r 
im 
1 
G 
r 
G 





Eur G. oF FG = ERGs) + EEGs 


wal PT AOS EE PETE O GT SRet lpi AG (4.2) 


V[FyGy + F4Gq] = E[(FyGy + FyGj)*] - E*[F,G, + F4Gq] 


PPE (Re {Cr FicG;- + 2FGFiGm) — 4pc6 dp> on" 


G 
See Er Ge Nererire CG; cl ttach le eGek. Gok Gp sO na Op 
= 9,7 spe ON Se Aa) + ee Cn lact 2p) 
4 20,7 akon a Po”) = Apo? a On? 
4 a2 gee? + 4p + Moe? - 2p" ~ 4p 7) 
= oar Ja a ee % Po’) (4.3) 


Similarly, for the imaginary part of the cross power 


spectrum, 


BG eG |i\= 2PQ On OG (AeA) 
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= L oe Nea 9 2 eae 12 
V[F.¢_ FG] BeBe el. Gs Pg one (4.5) 
Although the mean and variance of the spectra will be 
as derived above, the pdf will not necessarily be Gaussian. 
However, our concern is with the average of many 
instantaneous cross spectra, for which the pdf as shown 
below will tend to be Gaussian. 
Let the average over K instantaneous cross spectra be 
K 
2 


Px(m) = Po(m) + 5Pqim) = 2 1 F*(m,k) ee 


In the absence of interference, each instantaneous cross 
spectrum will be independent and by the Central Limit 
Theorem of statistics, P.(m) and Py (m) will tend to be 
normally distributed. The mean and variance of P, (m) will be 


K 
ee ieee ee eer ar ante 
1 § eae ae () 
Spyeicivecchateina C8 (4.7) 
ae 202(k) o2(k)[1 + p.2(k) + p,?(k)] 
VIP.(m] = yo 421 2%, G re PQ (4.8) 
Note that On» I » Po and Pg have now become functions 


of K in order to allow for the variation of these parameters 


Witthmesime. win practice, o and ». may change significantly 


C Q 
during one averaging interval determined by kK, due to the 
movement of discrete radio sources through the lobes of the 


antenna beam formed by the interferometer. OF and Oo wil] 
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change much more slowly, as they reflect the total power 
being received by the beam of each individual antenna from 
the entire sky. o, and So, may be considered to remain 
constant over the averaging interval. This is fortunate, as 
otherwise the requirement for identically distributed 
variables during summation for the Central Limit Theorem 
would be violated and the normality of the cross power 
spectra could not be expected. 

Taking °% and % as constant, 


K 


E[P.(m)} = 205, 6, = cape (k) = 25p % Pe (4.9) 
» ee : Z at 2(k)] 
VIP .(m)] = 20,7005? go 424 [1 tp? tk) + Pq (K) (4.10) 


Now Po and Pg were defined as the correlation 
coefficients of the instantaneous co and quadrature spectra. 


A complex correlation coefficient py may be defined as 
Py = Pot IPp (4.11) 


where the magnitude of Pe is equal to the fraction of the 
total power received by either antenna (assuming identical 
antenna patterns) which is correlated, O<|p, 151, and the 
phase Oe determines how the correlated power is 
distributed between the co and quadrature spectra. It 


follows that 


0 < py (k) py (k) ee cas Py (k) <a U4 412)) 


di? 


wae ae 
mot snneines ust VAP) 


wi 
-y 
A a } vw « Fe 
- ' a’ edt 
7 
ij } 7 
P a 
y! oe ! ty 
) 
rin hasieo s7 36 
fopue. SW RIOBUP OF OO aus 




















af re ro oe i) 


ting fato? at tot ar as ihe 
Pe 4 is. piate a 


; m4 
a di we 


a - 
erenut1a% et gin? . baicceth vie: 

ers 

-otuslinsath 41 Sok Mea AST Ral 


nfames of Be neti Me 


feesealh 
ot tangmae (pe 


me “Odi 7 an reine 2 ae fe 
1awoq gaot%D sft o ¥i tam _ “ ' 


25 beniteb ae vom 5 4 neta} tieea oe 


‘ : is 
) a ~ @ 


L e all ‘d - 
2 af oe : 


eit: Jo notiset? ets of Tape af ah” v3. sui tn = 


; niotek 
aaiineti pornimuses) she tite. peereye vs ; 


ons ee Poel Lottie PS. oo 


et young 'S sets sara tm 


| ii eh z enue ap be - 3 


ine: a 
aii, cs ra 7 a 


















o 
. n es ‘ 
»« “oe - or ae 
be >» ie oa ie eae ns 
pee. ©.) oe 
¥ 





vl a * 


er | oy 
gs " : 






= ta 


K 2 2 ; 
geen 2 [1 + poo Ck) + po? (k)] <a Av 3h 


and 20,° o,°/K < V[P(m)] s Wiice ep (4.14) 


At this point, an assumption may be made about 
decametric antennas which will place VIP. (m) ] close to the 
lower limit in equation 4.14. Practical decametric antennas 
generally have relatively broad beam widths, thus the total 
correlated power received from a discrete source wil] 
generally be small (lo, |<0.1) relative to background noise 
received from the rest of the sky. Hence spe (Ba Se Po USS 


and 


~ 2 
V[P (nm) ] fe 205, ore /K (4.15) 
The probability density function of P (m) is then found to 


be a normal distribution described as N(20, 0,0, oe an ei 


Similarly, the pdf of P.(m) is N(2 oo 459 » 


Q 20,7 9,7/K). 
4.2.2 Probability Density Function of Averaged Auto Spectra 
The instantaneous auto power spectra of the two signals 


will be 


FF*(m,k) = F(m,k) F*(m,k) = Fr? (m,k) fe F42(m,k) (4, 16) 


GG*(m,k) =G(m,k) G*(m,k) = Gr2(m,k) + Gy? (m,k) Ame 7a 


As a is independent of o and Ges independent of G. ron 
the case of no interference, the pdf's of FF’ and GG* will 
correspond to ae and OG respectively, times a chi-square 
distribution with two degrees of freedom, y?(2). The 
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averaged auto spectra will be 


P_(m) = ibis FF* (m i) (4.18) 
F K k=l : : 
G i K k=1 $ ; 19 


and will consist of the sums of the squares of 2K 
independent, normal variables. If o,? and o,* are assumed to 
remain constant over the averaging interval, then P.. (m) and 
P.(m) will have pdf’s corresponding to 0,?/K and 0, ?/K, 
respectively, times ¥2(2K), a chi-square distribution with 


2K degrees of freedom. For large K, »2(2K) will converge to 


a normal distribution N(2K, 4k), and so the pdf of P.(m) 





will be 
Op op. OR! 
sae NOTIN O) MOAN Freraty TAN rg oe aaa 4o,") (4520) 
K 
and the pdf of P .(m) will be 
o 2 
2 N(2K,4K) = N(20,2, 40") (And) 
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4.2.3 Estimating the Variance of the Cross Spectra 

In the discussion of robust estimation later in this 
chapter, it will be found that an accurate estimate of the 
variance of a distribution is invaluable when attempting to 
estimate the distribution’s center. It will now be shown 
that the auto spectra can be used to provide an accurate 


estimate of the variance of the cross spectra. 
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First of all, it can easily be demonstrated that the 
sample variance is a poor estimator of the cross spectrum 
variance. From Section 4.2.1, the pdf of P.(m) is 


N(20,0,0, » 20,7 g,?/K). The usual sample variance for P {m) 


would be 
= 2 
accu eamrce (4.22) 
aie M 
where Pes vl M) Dare Mies athe P(m) are independent for 


m=l1 C 
Gaussian noise, S? will have a pdf corresponding to 


of /(M-1) times v2(M-1), where g2 is the variance of Ben) 


The expected value of S* will then be o,? (because S? is an 


unbiased estimator of o,? [62]) and the expected variance of 


Sieewisltl abe 


lee 


Visite zor] ci] Ont 
a= rid (4.23) 


ba 


The expected percentage error in S? can be found by taking 
the ratio of the expected standard deviation of S? to the 


expected mean. 





2 
N[S?) 4 
% Error for Sample Variance = 100 - a = GNeo) 205, 1 
E{s*] a2 
M-1 C 
= 100 Z 
— (4.24) 


For example, if there are 100 spectral points, M=100, the 

expected error in the sample variance will be 

100( 2/99)=14.2%. Clearly, this error is quite large and a 

better estimate of the cross spectrum variance is desirable. 
A better estimate can be found by noting from equations 


Below e 4520) vand- 4321" that 
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a F2 , 
Tha 20, ats = (20,,°) (20,2) /2K 


B[P,(m)} + E[P.(m)]/2k = Pom) Po (m) /2K (4.25) 


The expected error of this estimate is determined by first 


finding the error in P,(m) and Py (m). 


PE (m) = P_(m) (4.26) 


As P,(m) is distributed as N(20,?, 40,4/K), the variance of 
P(m) will be (1/M)4o,4/K and the percentage ratio of the 


standard deviation to the mean will be 





4o_7 
ip F 
LaGrroreLor P_(m) = 100 [P.(m) ] = 100 MK 
F = a 
E[P_(m) ] 20-2 
F F 
100 (Am2/s 


A typical value of K such as 10,000 (corresponding to an 
averaging interval of 25 seconds) will result in a very 
small percentage error. The result for P .(m) is identical 
and the overall expected error in On*, which is derived from 
the product of two quantities with small percentage errors, 
will be twice that error or 


200 
mWuO bento .1 Cie = ae (4,28) 


For M=100 andk =10,000 the resulting error will be 0.2%. 
Even if a single point from each auto spectrum is used (i.e. 
M=1), the percent error is only 2%. Hence, the auto spectra 
can provide a much better estimate of the variance of the 


cross spectra than the sample variance. 
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4.2.4 Some Practical Considerations 

Two assumptions were made in Section 4.2.1 before 
deriving the normal distributions of the cross spectra and 
also the auto spectra. These assumptions must be justified. 

The first assumption was that the signals being 
correlated consist of Gaussian noise bandlimited to exactly 
the Nyquist frequency. Because ideal filters do not exist, 
in practice it is necessary to bandlimit an analog signal to 
less than one-half the Nyquist frequency in order to prevent 
aliasing. The spectrum of the signals is then not flat, but 
rolls off on the ends. In order to have a true normal 
distribution in the spectra, the frequency response of the 
receiving system and anti-aliasing filters must be flat. The 
simplest way to attain a flat spectrum is to discard the end 
points which are affected by the roll-off. 

The second assumption was that the correlation between 
the two signals is identical in amplitude and phase at all] 
frequencies. This requirement first of all means that the 
frequency response of the receiving system and anti-aliasing 
filters must be very flat in the center of the spectrum 
(disregarding the roll-off on the ends mentioned above), and 
that the phase response of the two channels must be very 
closely matched. Also, the propagation delays from the 
antennas to the correlator must be identical as a time 
difference will cause a linear phase shift across the 
spectrum. Finally, the signals received at the two antennas 


themselves must have a correlation which is identical in 
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amplitude and phase across the portion of the spectrum being 
observed. Identical amplitude of correlation is relatively 
assured because at decametric wavelengths the strength of 
Known radio sources changes very slowly as a function of 
frequency, and can be assumed to be constant over a smal] 
bandwidth such as 50 KHz. Identical phase, however, will not 
always be the case as there will generally be a time delay 
between the arrival of the signal at the two antennas. In 
Section 5.4 it will be shown that the maximum phase shift 
across the spectrum for the system during field trials was 
quite small and could be ignored. 

The measures taken to ensure the flatness of the 
spectra and phase matching are discussed in Chapter 3 on 
receiver design and Section 5.6 on field testing. These 
precautions assure that the cross and auto spectra will be 
normally distributed in the absence of interference, and 
allow the robust estimation techniques described in the 
remainder of this chapter to operate with maximum 


efficiency. 


4.3 A Survey of Robust Estimation 

4.3.1 A Brief History 

Perhaps since the time that the sample mean was first 
used as an estimate of the center of a set of measurements, 
practical people have realized that measurements which are 
grossly in error should not be included in the calculation 


of the mean. However, use of the sample mean has become so 
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entrenched that its users rarely consider that it is an 
optimum estimate only for errors which are independently, 
identically and normally distributed. Unfortunately, such 
normally distributed errors are often not the case. In the 
introduction to his survey of robust estimation, Ershov [63] 
mentions many reports of engineering measurements, 
industrial data, clinical medicine results and other 
situations where from 1 to 20% of data points may be 
considered to be anomalous, or "outliers" (i.e. not part of 
a normal distribution). In [64] Tukey quotes Geary: 
“Normality is a myth; there never has been, and never wil] 
be, a normal distribution." 

Gauss himself [65] introduced the normal or Gaussian 
distribution as that distribution of errors which was best 
suited in the sense of least squares to the sample mean, 
rather than vice versa. Before this, in the first published 
work on least squares, Legendre [66] notes that sample 
values which appear to be anomalous should be rejected 
before using least squares methods (and, in particular, the 
sample mean). 

A dogma of normality, and the sacredness of the sample 
mean, were helped by the Central Limit Theorem, which states 
that the sum of many small, independent elementary errors is 
approximately normal. Over the years, the sample mean has 
come to be used automatically in most situations, often 
without regard to whether or not it may be appropriate. 


Anscombe [67], on questioning the assumption of normality, 
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remarked: “The disposition of the present-day statistician 
theorists to suppose that all error Jeatr dbwenere are 
exactly normal can be ascribed to their ontological 
perception that normality is too good not to be true." 

Good statisticians and engineers have always been 
careful to exclude grossly outlying observation in mean 
calculations, and published accounts such as [68], [69], 
iO oma lll; elkf2'|, have proposed rules for rejecting outliers. 
But it was not until the 1960’s that a serious effort began 
to find estimation procedures with predictable properties 
which would provide immunity from erroneous observations. 

Section 4.3.2 outlines a history of rejection rules and 
literature on the performance of tests for outliers. Section 
4.3.4 describes work since 1960 in the area Known as robust 
estimation. Curiously, the literature on robust estimation 
does not generally encompass tests for outliers although the 
two topics are very closely related and overlap in some of 
their objectives. 

More detailed and quite interesting discussions of the 
historical development of ideas in robust statistical 
methods can be found in [63], [73], [74], [75], [76] and 
Ale 


4.3.2 Tests for Outliers 
In 1936, Pearson and Chandra Sekar [72] published an 
analysis of a rejection criterion attributed to Thompson 


[78]. Thompson's proposal was to reject observations Yee for. 
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where y and S are the sample mean and standard deviation and 
t ig chosen such that, for a normal population, the 
probability of rejection is small. This test has come to be 
Known as the Studentized deviate test. Pearson and Chandra 
Sekar point out a number of drawbacks of this test, namely 
that for small sample sizes it is only suitable if there is 
a single outlying observation and for larger sample sizes 
there is a fixed limit, depending upon sample size, to the 
number of outlying points which can possibly be rejected. 

The possibility is then suggested of applying the 
Studentized deviate test in a sequential fashion, as 
follows: (1) apply the rejection criterion to the n 
observations; (2) if k21 outliers are rejected, apply the 
criterion again to the remaining n-k observations, 
recalculating y and S each time; (3) repeat until no more 
observations are rejected. Pearson and Chandra Sekar note 
two disadvantages of this process: (a) if there are more 
outliers of a given amplitude than can be rejected for the 
sample size being considered, rejection will not occur and 
the process will stop; (b) care must be taken in choosing 
the rejection levels Lh ae etc., to reduce the risk of 
rejecting points which belong to a normal distribution. 

In 1950, Grubbs [79] proposed a test using the ratios 


of the sum of squares of deviations for a reduced sample to 
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and the Yc) are ordered sample values, Y(1) £Y(9) eh ey 
A similar test statistic D, is found by omitting the 
smallest sample Y (4) rather than the largest Yn) . Prescott 
[80] shows that Grubbs’ criterion is algebraically related 
to the Studentized deviate test and produces identical 
results, so in effect it is really the same test. 

Kudo [81] demonstrated in 1956 that the Studentized 
deviate test is an optimum test for the case where there is 
exactly one observation which deviates from a common normal 
distribution, the variance of which is unknown. Furthermore, 


if the variance o2 is Known, the optimum test is 
ly, - ¥[ > te (4.31) 


where +t is again a constant which depends upon the sample 
size and is chosen to keep the probablilty of rejecting 
points belonging to the normal distribution small. 

Other rejection procedures have been proposed but have 
been found to be unsatisfactory due to insensitivity or 
difficulties when more than one outlier is present. These 
tests include Dixon’s [82] criterion, David, Pearson and 
Hartley’s [83] criterion, and the skewness and kurtosis 


criteria of Ferguson [84]. A graphical comparison of the 
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behavior of all of the above tests when more than one 
outlier is present is given by Prescott [80]. 

The conclusion arrived at by Prescott is that the most 
satisfactory rejection prodedure for more than one outlier 
is a version of the sequentially applied Studentized deviate 
test introduced by Rosner [85]. For his method, one first of 
all selects a maximum number k of observations which one is 
prepared to regard as being outliers. The Studentized 
deviate test (or equivalently Grubbs’ test) is applied k 
times, each time rejecting the largest observation whether 
or not the rejection criterion is satisfied. The rejection 
procedure is then carried out in reverse, starting with the 
n-K remaining observations and reincluding one additional 
observation at a time, in the opposite order to that in 
which they were rejected, until an actual: outlier is 
identified. This point plus all of those beyond it in the 
order of reinclusion are then regarded as outliers. Rosner’s 
Methodwovercomes, satleast for-up to K outlying points, the 
Studentized deviate test’s disadvantage of not identifying 


any outliers if too many of the same amplitude are present. 


4.3.3 Rejection of Outliers in Robust Estimation 

Some statisticians working in robust estimation seem to 
harbor a strange distaste for rejection rules. Huber [73] 
cites a criticism of outlier rejection by Anscombe [86] and 
states: "The traditional philosophy behind rejection 


procedures is highly objectionable." Huber’s main objection 
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appears to be that the performance of rejection procedures 
is not easily amenable to analysis. Ershov, in his survey 
paper [63], says: "The fundamental deficiency in the 
procedures for rejecting anomalous data is the fact that 
they become complicated and not practically suitable for any 
significant amount of anomalous values (for example, 5 to 


20%), as well as distributions with heavy tails." However, 
Rosner’ s method (which Ershov includes as a reference) 
appears to be neither complicated nor unsuited to many 
anomalous values. 

Robust statisticians then proceed to include among 
their methods their own versions of rejection rules, under 
various guises, and to expound at length upon the merits of 
these procedures. A number of robust estimation techniques 


which may be called disguised rejection rules are described 


be low. 


4.3.3.1 The a-Trimmed Mean 

The a-trimmed mean consists of throwing away the [an] 
largest and smallest observations of a sample and taking the 
mean of the remaining observations. [| ] denotes the largest 
integer in the quantity in brackets, and QSas0.5. The 


estimate is then 


n-[an | 


1 : (4.32) 


Gre COs il) ake 
c n-2[an] 


SA) 


i=[an]+1 


where Y 4) are the ordered samples, Y1y S¥ (2) S++ SV (ay - 


The a-trimmed mean has a long history dating from its 
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use in certain provinces of France for calculating a mean 
annual harvest [87] to Poincaré [71]. Bickel [88] 
investigates the properties of this estimate and its close 
cousin, the Winsorized mean. The a-trimmed mean was included 
in the extensive study of robust estimation carried out at 
Princeton University in 1972 [89]. 

Major drawbacks of the oa-trimmed mean are that the 
trimming proportion @ cannot be chosen properly unless one 
has information on the distribution of contaminating points 
(outliers), and that because equal numbers of observations 
are trimmed from both sides of the sample the estimate does 
not perform well if the contamination is one-sided or 
asymmetric. 

If a has its maximum value of 0.5, the a-trimmed mean 
becomes the median of the sample. The median is a relatively 
robust estimate of the center of a contaminated distribution 
(it is much less sensitive to outliers than the mean). For 
this reason the median is often used as an initial estimate 
of location for the iterative robust methods to be described 


later. 


4.3.3.2 Skipped Estimates 

Tukey [90] proposed a class of estimates based on the 
rejection of outliers and called skipped estimates in the 
Princeton robustness study [89]. Some of these estimates 
were found to be among the best of the 68 different 


estimates considered but they have received suprisingly 
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little attention since that study. As an example of this 
type of estimate, the one designated 514 (one of the better 
ones) will be described below. 

First of all, from the ordered observations Yay? 
Y(1) SY (9) S.-Y (ny 5 the first and third quartiles are found. 


These are designated h, and h, - 


¥({n/4]) iedOteaasmulbtiplesok 4 
= (4.33) 
2 (a/4) na ¥ (14/4) » N a multiple of 4 
Y (nt+1-[ (n+3)/4)) » n not a multiple of 4 
Ss 
2 1 : : 
sh vette + Sreeainiy Area MULCrD Le vol 94 (4.34) 


hy is the value below which the smallest quarter of the 
observations lies, and h, is the value above which the 
largest quarter lies. The interquarterile range, h,-h, is 
used as an estimate of the scale (dispersion) of the data. 


Two points ba and t further from the center of the 


data than Me and h, are then found: 


x (Am35.) 
Peet eg coisa 


rt 
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a Rye Wen 486) 


The estimation procedure for 514 is as follows: (1) Delete 
all points less than t, and greater than t- (2, lt skipoints 
were deleted in (1), delete a further L points from each end 


of the sample where 


126 


Olek .— oe io - 



















git to ol quese ns 2A A vas ; 


ey 7! aot 8Y VWra5S boishten: ene a se . 
one doh ttamugl aha Abeta ext vi yes 
a gt tre b A toe 


rat Jar art? to 90} ATe hotel 


a 


. : a 
me 7 
\ if {car 6 Ug ih -» ; 
(ee 3 - 
; 7 irs N° « We \art 
; : b ‘ i 
| bhattiug & it 4 ( <l ‘oy 4 
(AER) ie \elarew) 
. he Fi i 
~ a4 - tr 4 * . a 7 = 
ie 3 | a a3 + 1405 A 29 ! tees 3! 3 Aol a woled il 
ont fotrew swoudr suf sv. ar ei A oils, EX = 
1 


at . fron .eenan pi + eitraduppacitred eat oa 9: 
-s36b sft ta toot paneqa hb} elsoe ‘ont *o"8 lemisi 


sti to isines: edie, nett? 7 oe sh 
Ornst nea one bane 
19 a ) y ar ip 
° 5 - ye 


7. 
. 
4 


ena! aes 
pean 9 ‘7 “i 


‘ ' ; a 7 +e nee Pact ; 7 
{ ae Bh) r : . € sh “jee 3 » 
: - - ; ais : 
— os . 
eraieal (it) : awed! hor ana a 

: Te . < 

sar W th (Slee tel sa 
tis se Lbee a A be 


¥ erie 25 We) s ee 

: ac . “ae i : i 
; ne, : JF) og. 
gare is “A 77 a 











L = min[max(1,2k), O23ne— k/2] AMES; 7s) 


(3) Take the mean of the remaining samples as the estimate. 

Though it is very simple and easy to calculate, even by 
hand, this estimate was found to perform well in the 
Princeton study and later in a study by Wegman and Carrol] 
[91]. It is superior to the a-trimmed mean because the 
decision of how many points to delete is made on the basis 
of the number of outlying points detected in the sample 
itself, rather than being fixed beforehand. The 
disadvantages of the 57T4 estimate are that the interquartile 
range, used as an estimate of scale, is not the best (see 
Hampel [92]), and also that the estimate of the center of 
the distribution used in rejecting outliers, namely the 
center of the interquartile range, is also not particularly 
good. In addition, the rejection of an additional L points 
based on the number kK originally deleted may or may not be 
appropriate. It is not surprising, then, that in [91] the 
51T4 estimate was generally inferior to another class of 


estimates called M-estimates, described in Section 4.3.5. 


4.3.4 Maximum Likelihood Est imation 

One of the most important estimation techniques in 
statistics is the method of maximum likelihood reviewed in 
this section. 


The likelihood L of an observed sample nots cass Canc a 
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from random variables Ky X ae x is defined to be the 
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random variables. A maximum likelihood estimator produces an 
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estimate of a parameter such that the likelihood (i.e. joint 
probability or joint density) of the sample is maximized. 

To illustrate, assume some density function f(x) is 
displaced by an unknown amount oe to form a density f(x-e). 
To find a maximum likelihood estimate of the location 
parameter 6 from a sample of independent observations Xo 


x x0 first of all note that the likelihood function 


pe et ae 


L(e) will be the joint probability of x,, x for 


; ey eX 
Hh 2 n 
some value of 6, and will be given by the product of the 


individual probabilities. 


L(@) = ,T, £(x,-8) (4.38) 


Because it is simpler to deal with sums rather than 
products, it is often convenient to maximize the logarithm 


Ofuslcuerathner ethane |.o) ma tselt. 


n 
a 


Ins C6) ; 


1 In f(x ,-6) 
n 


= -,E, p(x,-8) (4.39) 


where p(x)=-In f(x). The maximum is found by setting the 


derivative with respect to 6 equal to 0 and solving for 6. 
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where w(x)=-p’ (x). The solution of 


nh 
2, p(x,-@) = 0 (4.41) 


1 


which maximizes L(6) is the maximum likelihood estimate of 6 
and is frequently denoted as 6. 

It is easily shown that the maximum likelihood estimate 
for the normal distribution is the sample mean. Dropping 


some irrelevant multiplicative constants gives 
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x2 
pO) = 5 
Wx) = x 
n n 
2, WOs,-8) = 8B, (x,-6) = JE) x, - n@ = 0 
re | n — 
Therefore, 6i= = way us was 


The importance of the function w(x) for robust maximum 
likelihood estimation will become clear later. For now, 
consider the significance of p(x). First of all, the sample 
mean is Known as a least squares estimator because p(x) has 
the form o(x)=x2/2. The estimate 6 minimizes the sum of the 
squares of the deviations from 6, i.e. 6 maximizes 


n 


52, p(x,-8). A maximum likelihood estimator may also be 


thought of as maximizing a correlation between the 
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non-displaced sample density function f(x) and the function 
-p(x). 

To illustrate this idea of a correlation being 
maximized, consider the case of an infinite sample size, 


n+”, Maximizing the log of the likelihood function, 
nN 
In L(6) = -,2, p(x,-8) (4.42) 


as n+- will be accomplished if the expected value of 
OX ahr.ot) is maximized. xX, may be replaced by a continuous 


variable x, and the expected value is 
E[-p(x-8)] = J -p(x-e) f(x) dx (4.43) 


The integral has the form of a correlation, R(e), between 
the function -p(x)=In f(x) and f(x). The correlation will 
have its maximum value for 6=6. 

The purpose of demonstrating this connection between 
maximum likelihood estimation and correlation is to contrast 
the estimation problem to the optimum receiver problem in 
communications theory. It is well Known [93] that the 
optimum receiver for a signal in the presence of noise is 
one which correlates the signal plus noise with a duplicate 
of the signal alone. The correlation will then be a maximum 
compared to the correlation with all other signals of equal 
energy. However, for maximum 1ikelihood estimation one 
attempts to maximize the correlation between a probability 
density function f(x) and its logarithm rather than a copy 


of f(x) itself. For example, if one correlates a normal 
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distribution with a -e(x) having the same Gaussian form, the 
resulting estimator performs rather poorly (c.f. the MEL 
estimate used in the Princeton study [89]). There is an 
intrinsic difference in the nature of the two problems. The 
optimum receiver problem attempts to detect the presence or 
absence of a particular waveform, whereas the estimation 
problem attempts to determine the location of a density 


function as accurately as possible. 


4.3.5 Robust M-Estimators 

One of the major advances in robust estimation occurred 
in 1964 when Huber published a paper [94] introducing a 
class of maximum likelihood type estimators, subsequent ly 
known as M-estimators. These estimators are solutions 6 to 


the general maximum likelihood estimation equation 4.41, 
n 
i271 v(x, -6) =a) 


for various forms of the function wv. Huber showed that under 
quite general conditions (normally encountered in practice) 
f converges to an asymptotic value 6, defined by 


co 


CXS ORE Go) x= 20) (4,44) 


- oc 


for a given density function f(x). If f(x) is symmetric and 
w(x) is an odd function, then @, is the center of filaxole 
Also, as n>, Gioseal ce 8 ay is asymptotically normal with 


asymptotic mean 0 and asymptotic variance 0 given by 
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62 ne [p(x-6_) 1° £(x) dx 
(4.45) 


co 


[s, W (x-0,) f(x) dx}? 


Scale invariant versions of M-estimators are obtained 
by dividing by the standard deviation 9° (or an estimate of 


c). The estimation equation then becomes 


o ss] 
a vf |- 0 (4.46) 


O 





This equation is usually solved iteratively from some 


initial estimate ®,of ® (such as the median) using the 


Newton-Raphson method. Thus, 








(3) |e 
k+l k A x ,-8 
S| seLwak (4.47 


Some of the commonly used ~ functions recommended by 


Hogg [95] are given below: 


1. Huber Proposal 2 
W(x) =(-k, x < k 


Say ties (4.48) 


A good value for kK is 1.5. 


2. Hampel 
Ix] , OP <aaix meuta 
W(x) = (sign x). a, B4< ope ler< ib 
a C—7 |X| CD < mice ice 
c-b 
0: Ors) 1x (4,49) 
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Reasonbly good values of the constants are a=1.7, b=3.0, 
and c=6.4 which correspond to an estimate Known as the 
Hampel 25A. 


3. Wave of Andrews 


sin (2%), ball xe ik 
£ k 
fey e 0 ont oak (4.50) 


with k=5.0 or 6.0. 


4. Biweight of Tukey 


2 
Voy eye ata , [xl sk (A541) 


0 exe ik 


with k=5.0 or 6.0. 

Because the last three w functions above go to zero for 
|x| beyond specified limits (a desirable property because 
extreme outliers are rejected completely), there may 
sometimes be convergence problems with their associated 
estimates. For instance, if the initial estimate of 6 used 
is not close to the bulk of the data there may be too few 
points within the limits of the » function for proper 
convergence. These problems can generally be avoided by 
assuring that the initial estimate of 6 is reasonably close 
to the final value. Hogg [95] recommends that Huber’ s 
function (which does not suffer from convergence problems } 
with k=1.5 be used for several iterations to provide an 
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4.3.6 Other Robust Estimators 

Besides M-estimators, there are two other general 
classes of robust estimators called L-estimators and 
R-estimators. L-estimators consist of linear combinations of 
order statistics (i.e. observations are weighted according 
to their rank rather than their actual value). R-estimators 
are based entirely upon the ranks of observations rather 
than their actual values. Both of these classes of estimates 
are described in [89], [73], [63] and [95]. It has been 
shown ([96], [97]) that M-estimates can generally be found 
which are asymptotically equivalent to L- and R- estimates. 
Also, in Monte Carlo studies such as [89] and [91], 
M-estimates are generally found to be superior to the L- and 
R-estimates studied. Therefore, these two classes of 


estimators will not be considered further herein. 


4.4 The Performance of M-Estimators 

4.4.1 The Influence Curve 

An important tool for studying robust estimators is the 
influence curve developed by Hampel [92]. The influence 
curve is essentially the first derivative of an estimator 
evaluated at some distribution, and can be used to derive 
asymptotic variances and several other robustness properties 
which will be explained below. 

To define the influence curve, let F(x) denote a 
probability distribution and T(F) denote some estimation 


function defined for a subset of all probability 
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distributions, including F. Consider a probability 
distribution consisting of a mixture of F and a delta 
function at some point x, written as (1-c)Fted,, O<e<1. Then 
the influence curve IC, (Xx) of the estimator T at the 


under lying distribution F is defined pointwise by 


TC) p(x) = lim ii (l-e Fre é, |=) TC) 
€>0 € (4a5 2) 


The value of the influence curve at a point x is equal to 
the change produced in the estimate T by the addition of a 
point mass 1 of value x to the underlying probability 
distribution. 

For M-estimators, the influence curve is very closely 
related to the w function. In particular ([92], [89]), if 


f(x) is symmetric and ¥ is odd, then 


Bee) feet e (4.53) 


where the term in square brackets is just a scaling factor. 
The asymptotic variance in this case can be found from the 
expected value of the square of the influence curve. 

02 SF [w(x/o) ]2£ (x) dx 


Ceegr (LC a) ge = (4.54) 
[Jw (x/0) f (x) dx]? 





Influence curves for the four M-estimators from Section 
4.3.5 are shown in Figure 4.1. A fifth M-estimator, 
corresponding to a rejection rule which deletes points 


beyond +ko (k=3.0 in the figure), is also shown for 
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FIGURE 4.1 Influence Curves for Robust M- Estimates 
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4.4.2 Properties of M-Estimators 

Three important properties of robust estimators can be 
derived from their influence curves. These properties are 
gross-error-sensitivity, local-shift-sensitivity, and 
rejection point. 

Gross-error-sensitivity measures the worst possible 
effect which a single contaminating point can have on the 
value of an estimate. It corresponds to the maximum absolute 


value of the influence curve, and will be denoted GES. 
GES = max eese (4.55) 


Good GES is obtained by placing a bound on the influence 
curve and making that bound as small as possible (the 
influence curve of the mean, in contrast, is IC(x)=x and is 
unbounded). In general, the goal of limiting the influence 
curve (and hence limiting errors due to outliers) conflicts 
with the goal of having an asymptotically efficient (i.e. 
smallest variance) estimator. As the bound is made smaller, 
the efficiency relative to a most efficient estimate such as 
a maximum likelihood estimate generally decreases. 
Fortunately, the price which must be paid in efficiency is 
very small (a few percent) to obtain far greater gains in 
protection from the effects of outliers. The beneficial 
trade-off is what makes robust estimators practical and 


attractive. 
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A second property of robust estimates is called 
local-shift-sensitivity, denoted as LSS, and is a measure of 
the worst possible effect on the estimate of a small change 
in the value of a single observation divided by the size of 
the change. The LSS is found to be the maximum absolute 


value of the first derivative of the influence curve. 
LSS = max i IC(x) | (456) 


A large value for LSS does not imply that an estimate wil] 
necessarily be inaccurate, but means that the estimate is 
very sensitive to changes in the distribution of 
observations near some particular values of X and as a 
result could behave erratically. This fact is particularly 
relevant when one considers that for a finite scale size 
the distribution function will not be smooth but can exhibit 
local grouping of points. Contamination may also cause loca] 
grouping. A very high local shift sensitivity turns out to 
be the major disadvantage of rejection of outlier estimates. 
The third property of estimators is the rejection 
point, which is the point beyond which the value of the 
influence curve becomes zero. All observations beyond the 
rejection point are rejected completely. A low rejection 
point is desirable to eliminate as many outliers as 
possible, but its attainment conflicts with the requirements 
for high asymptotic efficiency and a smal | 
local-shift-sensitivity. Generally, a compromise must be 


reached. 
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4.4.3 Comparison of Some Estimators 

The numerical properties of the M-estimates discussed 
above for a standard normal distribution are summarized in 
Table 4.2. The mean and median are included for comparison 
purposes. Figures for estimators a, b, d, and e were taken 
from [92]. No published figures for the remainder were 
available, so they were calculated numerically by computer 
using equation 4.54 for o, ? and equation 4.53 for the 


influence curve. 


4.4.3.1 30 Rejection Rule 

A 30 rejection rule, which takes the mean of all 
observations within +30 of the center, has low asymptotic 
variance and a low rejection point. However, GES is the 
largest for any of the estimates considered and LSS iso, 
indicating extremely high sensitivity of the estimate to 
observations close to the rejection point. This estimate 
could be expected to behave badly in the presence of 
low-level contamination, but would perform well for 


contamination greater than 3c. 


4.4.3.2 Huber’s M-Estimator 

The Huber estimate has low asymptotic variance, low 
GES, and low LSS. In fact, the Huber estimate can be shown 
to have the smallest asymptotic variance of any estimate 
with a given gross-error-sensitivity (i.e. it is a minimax 


estimate [94]). Its disadvantage is an infinite rejection 
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Table 4.2. Numerical Properties of Some Robust M-Estimates 


Estimate 
a) Mean 
b) Median 


c) 30 Rejection 
d) Huber k=1.5 
e) Hampel 25A 
f) Wave k=5 

g) Biweight k=5 
h) Biweight k=4 
i) Biweight k=3 
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point, which as mentioned in Section 4.3.5 ANS vie ee 


convergence to the true center of a distribution regardless 


of the initial estimate used, but which degrades its 


performance in the presence of large contamination. 


4.4.3.3 Hampel’s Redescending M-Est imate 


As seen from their 


is simply a Huber estimate which redescends to zero rather 


than staying constant. 


inf luence curves, 


a Hampel estimate 


Thus, a Hampel estimate overcomes the 
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Huber’s deficiency by having a finite rejection point (at 
the expense of slightly increased oh for the same GES). 
Hampel estimates are therefore recognized as being good in 
all respects ([89], [92]). By varying the constants a, b, 
and c various performance factors may be traded off with one 
another. For instance, the rejection point may be improved 


at the expense of some of the other factors. 


4.4.3.4 Wave M-Est imate 

The use of a smoothly changing wv function such as the 
biweight or wave has some slight advantages over one 
composed of straight line segments like the Hampel. For 
instance in Table 4.2 the wave has the same variance and a 
slightly higher GES than the Hampel, but it must be noted 
that the region for which the influence curve is a maximum 
is much smaller than that for the Hampel. The wave has a 
lower rejection point which is paid for mainly by a larger 


ESS. 


4.4.3.5 Biweight M-Est imate 

The wave and biweight estimates are very similar and as 
Hogg [95] suggests, almost interchangeable. The performance 
of the biweight with k=5.0 is very close ot that of the wave 
with k=5.0. Both of these estimates have very good qualities 
and in the opinion of the author of this thesis are the best 


choices for most applications. 
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4.5 Adaptive Estimation 

Three versions of the biweight are shown in Table 4.2 
with values of the parameter k of 5.0, 4.0, and 3.0. Note 
that as the rejection point is decreased, the variance cha 27 
increases. LSS also increases but GES is actually reduced 
somewhat . 

The asymptotic variance oe is the variance for very 
large sample sizes (n+~). Because the actual variance of any 
estimate varies with the sample size as 1/n, for n>= the 
actual variance should approach zero. The interpretation of 
one in Table 4.2 is that shite is the variance for a very 
large but finite sample size, relative to the variance of 
the mean (the optimum estimate for a normal distribution). 
If one is dealing with a large enough sample size the 
expected variance may be small enough that one would, in 
some situations, tolerate an increased variance in order to 
gain greater protection from outliers which are relatively 
close to the center of the distribution (within a fewo). 
Such a situation would occur if there was evidence or 
suspicion of low-level contamination, and one wished to 
minimize its effects on the estimate. In such a case, a 
biweight parameter of k=4.0 or 3.0 might be more appropriate 
than k=5.0. 

Adaptive estimation consists of choosing which estimate 
to use after looking at the data and making some decision 
about the nature of the contamination. If one could through 


some choice of reliable statistics determine that there was 
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a significant amount of close-in contamination and could 
then choose the value of biweight parameter which would 
produce the most accurate estimate in the presence of this 
contamination, then adaptive robust estimation would be a 
major improvement over estimation based on a single biweight 
» function. 

A few relatively simple forms of adaptive estimators 
have been suggested and evaluated in [89], [98], [99], [91], 
and [100]. These estimators have not been spectacularly 
successful, for though they may be better than non-adaptive 
estimators in worst-case situations, they are generally 
somewhat worse in "run-of-the-mill" cases. The problem 
appears to be a difficulty in judging the nature of close-in 
contamination. Often the adaptive prodedure will select a 
non-optimum y»-function (i.e. one for more contamination than 
is actually present), and consequently overall performance 
tends to suffer. However, most authors express much hope 
that adaptive prodedures can be improved (little effort 
seems to have been expended in this area as of yet) and so 
adaptive estimation should bear considerable attention in 
the future. Adaptive estimation would be of the greatest 
benefit in cases of asymmetric contamination. If outliers 
are concentrated on one side of the center of a distribution 
but are close enough to not be rejected, they will cause the 
estimate of the center to be biased in the direction of the 
contamination. A natural criterion for judging the 


performance of an estimate is the mean squared error 
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considered by Jaeckel [96], who shows that the mean squared 
error with asymmetric contamination can be minimized by 
trading off increased variance for decreased bias. This 
could be accomplished by reducing the value of the parameter 
kK for a biweight estimate. 

Asymmetric contamination is a major problem in the 
developing theory of robust estimation. An observation by 
Hogg is mentioned in [91]: "In particular, as Hogg has 
suggested in private communication, although a sample may be 
drawn from a symmetric population the sample may have 


significant asymmetries." For the rejection of terrestrial 
interference considered herein, asymmetry is very likely 
because one contaminating signal may affect a large number 
of points in the spectrum, first of all through leakage and 
secondly through a large occupied bandwidth (e.g. amplitude 
modulation). As of yet, no consensus has been reached on how 
robust estimation should handle asymmetric situations. Some 


discussion of the problem is contained in [96], [101] 


ate SUIRORL ED 


4.6 An Estimation Procedure 

There is no one estimation procedure which is optimum 
POpecleeciMi di} ONS ebacieOr, the methods ‘described in Section 
4.3.5 has some disadvantages. A combination of procedures, 
OamyOmid saporoach. can be used to overcome some of the 
individual disadvantages and is most appropriate where a 


wide range of contaminating signals must be handled. 
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However, even a very complex hybrid procedure is unlikely to 
work in all conceivable cases. A practical compromise is to 
keep the estimation procedure as simple as possible, yet 
able to handle the most commonly encountered contamination 
in a highly efficient manner and less common contamination 
with reasonable efficiency. Provisions should be made for 
difficult cases to be at least recognized and perhaps 
recorded. In this way one could evaluate the estimator’s 
performance and note improvements which could be made in the 
handling of difficult contamination. 

The estimation procedure described below is not claimed 
to be the best possible, nor is it the only way to achieve 
comparable results. Rather, it is one of many possible 
combinations, but one which the author believes will work 
well for the problem of contamination of the radio spectrum 
with terrestrial interference. Possible additions and 
improvements to the procedure are mentioned in Section 


4.6.5. 


4.6.1 Estimation of Variance 

As described in Section 4.2.3 the variance of the cross 
spectra may be estimated from the auto spectra. The auto 
spectra will be contaminated by the same interference as the 
cross spectra, hence some form of robust estimation is again 
required. The first decile (see below) is proposed as a 
simple but effective estimator. 


The auto spectra are different from the cross spectra 
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in that the auto spectra are always positive and the 
contamination is entirely one-sided. Interference may 
increase the power received at some frequencies, but can 
never reduce the power below the level of background noise. 
With no interference, the dispersion of the auto spectra 
will be very small relative to the mean. From equation 4.20 
the probability density with no interference (and assuming a 
perfectly flat frequency response) will be N(2o,?, foe? [Ky 


The ratio of the standard deviation to the mean is then 


op 4a TR 


y 


Hy 


= 1 


Ree 3 FE (4.57) 


a 


which will be .01 or. less for a typical kK of 10% or more. 
Almost any estimate of the mean, as long as it is not 
drastically affected by interference, will provide 
reasonable accuracy. 

In addition, the variance of the cross spectra need not 
be Known with extreme accuracy. Good robust estimators 
(especially those with small local-shift-sensitivity) are 
relatively insensitive to small changes in the estimate of 
scale used (this corresponds to the notion of "qualitative 
robustness" discussed in Hampel [92]). Accuracy of the scale 
estimate to within a few percent would appear to be more 
than adequate. 

Te estimator proposed for the centers of the auto 
spectra is the first decile (i.e. the point below which 10% 


and above which 90% of the observations are found). This 
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estimate is robust in the presence of single-sided 
contamination and, as will be shown, has sufficient accuracy 


for the required purpose. 


4.6.2 Accuracy of the First Decile as an Estimator 

Let A represent some fraction, 0<d<1, of the cumulative 
area under a standard normal curve, f(x). Let a, represent 
the value of x at which the area has been accumulated, as 


shown in Figure 4.2. 


oe 
A = im ‘ba exp(-x2/2) dx (4.58) 
If n independent observations, Yar Voreree YQ from an 


N(0,1) distribution are ordered to form n order statistics, 
Y 1) $¥ (2) S$: SV py then the expected value of the k’ th 
order statistic,y,, | will be a, as above with meio hy 2h talc 
It is Known [103] that the asymptotic distribution of order 
statistics approaches a normal distribution given by Na, 
MploxeeOT 24a ))). The mean and standard deviation of this 
distribution is evaluated for n=100 and various values of 4 
in Table 4.3. 

When using an order statistic P ik) corresponding to the 
first decile (k=[An]+1) to estimate the center of the auto 
spectra, the estimate will from Table 4.3 typically 


underestimate the mean by 1.3 oa and have a standard 


P 
F 
deviation of 0.17 o, +» From equation 4.57, the ratio of the 
F = 
standard deviation o, to the mean P, is 1/V7K . A small 


F 
correction factor of 1/(1-1.3/K) will correct for the 
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FIGURE 4.2 Cumulative Area Under a Standard Normal Curve 
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Table 4.3. Asymptotic Expected Value and Standard Deviation 
of Quantiles of the Normal Curve 


Quantile Expected Standard 
Value Deviation 
r on £-1(a,) YA (1-2) /n 

Ona eis! Onele7 

0.2 NOSE tebe) 0.14 

O25 0.0 Gh shee) 
underestimation and provide an unbiased estimate P.__ of the 
mean P. Thus 

paws 1 P 
EST ——— +1) 
(1-1.3/VK) Sot 
= (14+1.3/VK) Prandtl) fOreke ot (4.59) 


The standard deviation ofaee on will be 


Oeil tbg3/VK) $0217 fe = (14+1.3/¥K) 0.17 P/¥R 


= 0.0017 P for K = 10% (4.60) 


If no interference is present the first decile plus the 
correction factor above provides a very accurate estimate of 
the means of the auto spectra. 

If there is interference, the estimate is affected only 
Slightly. For example, if 50% of the points of an auto 


spectrum contain interference and are larger than they 









naiteived Gawbra te Yew soe 


avd a 


tn t2 a” 
au eye ee: 7) 
Sf Ved 
a j rt. yy" 7 
i 
r 
~ 6 
Gow: i! 











“2 @} re 
| Yineees 
wi . “+ nl fie a tong Lite 


ANE SaTy SS 


q ane. 10) Bo . 


—_ Me 7 


og Pree, ia iD nord ae 5 DIBA. 


+i nal) 


tr bes 
‘ 
t ie ans y Wit i+ ‘ a 4 ade 4 
L ; ; : a 
a 


re AO 4 te 48%, y 


ath ae ay a 






Le : = is 


ont eu be af 3 toes Jee 


piers 
a 


ip 
bay 
“in 


air ise ST eNUESE. eee a 
“ : ‘J/ 


should be, then the remaining 50% will still belong to a 
normal distribution and the order statistic + dia kab which 
previously corresponded to the first decile will now 
correspond to the second decile (\=0.2) of a normal 
distribution. The new estimate with the old correction 


factor will be 


P = (151. 3/ VK 
pop) ttt. 3/¥K) CTOs 28D 


Cal WES) (GIS SIE. ES than (4.61) 


for K=104. The error with 50% contamination is only 0.44%. 
The first decile is therefore insensitive to even large 
amounts of interference and will serve as a very simple but 
accurate estimation for the means of the auto spectra, and 
hence the variance of the cross spectra. 

It must be noted that this estimator will be affected 
by any anomalies, especially dips, in the frequency response 
of the auto spectra. The spectra must either be flat or any 
ripples must be corrected for before the first decile will 
function properly as the auto mean estimator. 

An alternative to the first decile which would provide 
better accuracy would be an M-estimator such as a biweight. 
However, the additional computational complexity of an 


M-estimator is not warranted in this case. 


4.6.3 Location Estimation 
The estimation procedure used to find the location 


parameters (centers) of the cross spectra described below is 
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a hybrid of three different procedures. A Huber estimator 
with k=1.5 is combined with a 60 rejection rule to obtain a 
reliable initial estimate for a final biweight M-estimate. 
The rejection rule takes advantage of expected sidelobe 
levels to reject additional contaminated points on each 
iteration, and combines rejections from both the co and 
Quadrature cross spectra. 

Prior to estimation, all spectra must be as flat as 
possible. First of all, end points affected by roll-off are 
deleted. Ten points are deleted from each end leaving a 
total of n=108 spectral components of the original 128. The 
auto spectra are corrected for amplitude ripples, and the 
cross spectra are corrected for amplitude ripples and, if 
possible, phase ripples. 


The steps in the estimation procedure are as follows: 


1. Order the components of the auto spectra according to 
magnitude. Find the location parameters of the auto 
spectra, o.. and 80) using the first decile plus the 
appropriate correction factor as the estimate. The order 
statistic which most closely approximates the first 
decile is found using 


k=[An]+1=[(0.1) (108) ]+1=11 


where [ ] is the largest integer function. Thus, 


Ch = Ie Gea (4.62) 
GL) 
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The standard deviations of the cross spectra are then 


estimated as 


o 
Ce Ke (4.64) 


Find the medians Mo and Us of both the co and quadrature 
spectra. Do not disturb the order of the spectral 
components in doing so. Use the medians as initial 
estimates of the location parameters of the spectra 


mime ale t .08=M ee andicd =MQ)- Calculate the standard 


Cig2G Q 
deviations about the location parameters: 
nN 
AR Ses YORMEET pats ee 
cme Had Miata (4.65) 
nh 
=. (2. (4 a1 1/2 
5 = [si g(-2)" (4.66) 
Q n 


If both S($1.25¢, and Sg81-250,, conclude that there is 
no significant interference and proceed to biweight 
estimation, step 10. 
lf either So or 5 is greater than 1.250,, outliers may 
be present. Starting with the spectrum (co or 
quadrature) with the largest standard deviation, proceed 
to delete outliers. Find the spectral component 
corresponding to the largest standardized deviate from 
the location estimate: 

Ge = max IPG) -e,| 


EES TS TS EN (4.67) 
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IF D ax to 9» delete this observation from the sample. 
Delete adjacent spectral components which will also be 
contaminated due to leakage. If Des 210.0, delete two 
adjacent observations on each side. Otherwise, delete 
one adjacent observation on each side. 

Repeat step 3 a maximum of 5 times on the reduced 
sample. Stop after 5 iterations or if no more points are 
found with D26.0. 

For the other spectrum, first of all delete the spectral 
components corresponding to those deleted in steps 3 and 
4 from the first spectrum. Then proceed as in steps 3 
and 4 to delete additional points. The combined 
components rejected from both spectra upon completion of 
this step are considered to be contaminated with 
interference and are subsequently rejected from both the 
co and quadrature spectra. | 

Use a Huber M-estimator (equation 4.48) with k=1.5 to 
estimate the location parameters of the remaining 
components in both spectra. Use the Newton-Raphson 


method of 4.47. 
: n _ [Pc(i)- a] 
Te a a S,, oR vES fbx 4 
3 y aa ex (4.68) 
i=l a 
Exclude deleted observations from the summations and 
start with the medians from step 2 as initial estimates 
of 8. Iterate until either the change in ® is less than 


some fraction of O,. (.01, for example) or until some 


fixed number of iterations have been completed. The 
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resulting estimates of 8. and 8 should be reasonably 
accurate. 

If either of the estimates has changed by more than 

0. lo, , indicating considerable asymmetric contamination, 
go back to step 3 and repeat the rejection and Huber 
estimation using the new estimates of Fe and 8Q° 
Reinclude all observations prior to starting rejection, 
as some previously rejected points may now not be 
identified as outliers. Continue steps 3 to 7 until the 
estimates become stable. This should normally happen 
after very few iterations. 

Calculate the standard deviations of the remaining 
points after deletion as in equations 4.65 and 4.66. If 


both S($1.250, and S,$1.250,, conclude that no 


Q 
Significant interference remains undeleted and proceed 


to biweight estimation, step 10. 
If either So Ores 


Q 
interference which should be deleted may remain. Note 


are larger than 1.256, some 


whether or not the maximum of 5 outliers was detected in 


steps 4 and 5 for either the co or quadrature spectrum. 
If less than 5 were detected, proceed to step 10. 
Otherwise, return to step 3 and proceed to delete up to 
another 5 points (plus adjacent observations) from both 
spectra. Continue on from step 4 to recalculate new 
estimates of 6, and ®.. Repeat steps 4 through 9 until 


C Q 
either (a) the criteria in step 8 are satisfied or (b) 


no more points meeting the rejection criterion in step 3 
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remain. 

10. Use a biweight M-estimator with k=5.0 to find final 
estimates of 8. and 8 from the spectral components 
remaining after the largest outliers have been rejected 
in steps 1 through 9. Starting with the last Huber 


estimate from step 6, use Newton-Raphson iterations 


until ® changes by less than 0.016,. 


4.6.4 Evaluation of the Estimation Procedure 

The above procedure is quite conservative and wil] 
operate successfully to minimize the effects of interference 
in the majority of cases. Steps 1 to 9 together provide a 
very good initial estimate for step 10, a biweight 
M-estimate. A Huber estimator is used in step 6 because of 
its good convergence properties, while the iterative 
rejection procedure removes large and easily distinguishable 
outliers to overcome the Huber estimate’s infinite rejection 
point. 

The procedure should not break down except in extreme 
situations. One such possibility would be one-sided 
contamination of more than 50% of the spectral components, 
in which case the median would be unreliable as an initial 
estimate of location. With symmetric contamination the 
procedure should work properly with even more than 50% of 


the points being contaminated. 
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4.6.5 Possible Improvements 

Two possible improvements would be to provide a more 
reliable initial estimate than the median to reduce the 
likelihood of breakdown when more than 50% of points are 
contaminated, and to use some form of adaptive estimation to 
reduce the mean squared error when considerable low-level 
interference is present. 

One way of increasing the reliability of the initial 
estimate would be to note that there is a certain range 
beyond which the centers of the cross spectra are not 
expected to be found. Other a priori information, such as 
the values of estimates for previous spectra, might also be 
used. A correlation of the histogram of the cross spectra 
with the expected normal distribution could be employed to 


initially locate the center of the distribution. 


4.7 Summary 

The cross and auto spectra are shown to have normal 
probability distributions with interference causing 
contamination and therefore heavier tails. The expected 
means and variances of the underlying normal distributions 
are found. It is shown that the location parameters of the 
auto spectra will provide a reliable estimate of the 
variance of the cross spectra. 

Published literature in the areas of rejection rules 
and robust estimation is reviewed. M-estimators are singled 


out as the best class of estimators, and properties of a 
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number of M-estimators are described. 

The order statistic corresponding to the first decile 
is proposed as a simple but effective estimator for the 
location parameters of the auto spectra and thereby the 
variance of the cross spectra. The performance of this 
estimate is discussed. 

An estimation procedure for the centers of the cross 
spectra consisting of a combination of an outlier rejection 
test, a Huber M-estimator and a biweight M-estimator is 
presented and evaluated. Some possible future improvements 
to this procedure, including the addition of an adaptive 


biweight M-estimate, are suggested. 
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5. Field Trials 


5.1 Introduction 

Field testing of the complete interference excising 
system was carried out between Nov. 12, 1979 and Jan. 6, 
1980 at the Dominion Radio Astrophysical Observatory (DRAOQ) 
in Penticton, British Columbia. Reasons for the choice of 
DRAO were the availability of a large 22.25 MHz T-array 
radio telescope and the expertise and experience of DRAO’s 
staff from a long-standing program of research in 
low frequency radio astronomy. 

The timing of the observations, by chance rather than 
design, coincided with an approximate maximum in the 11-year 
sunspot cycle. Due to increased solar and hence ionospheric 
activity, terrestrial interference and also scintillation 
and refraction of radio signals were expected to be quite 
bad. 

The part of the year from November to January is ideal 
for low frequency astronomy as a result of the long winter 
nights. Electron density in the F layer of the ionosphere 
reaches a low value a few hours after sunset and remains low 
until sunrise. Hence there is a long period during the night 
with low levels of terrestrial interference. During the day, 
however, electron densities reach maximum values (due to the 
winter anomaly) and so daytime interference is at its worst. 
The diurnal variation in electron density, with steep 


gradients at sunrise and sunset, is reflected in graph (d) 
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Ofehrounem ia. I. 

The interference excising system was not expected to 
handle the very high levels of interference during the 
winter day, and so the periods of transition at sunset and 
sunrise were of particular interest as a test of the system. 
Late fall was deemed a good time for testing because the 
strong radio source Cassiopeia A has an upper transit in the 
evening and a lower transit in the morning, and could be 
expected to provide easily observable fringes from which 


interference could be excised at these times. 


5.2 Antennas 

The antennas were originally DanteOtfedec2..25.NHz 
T-shaped array telescope at the Dominion Radio Astrophysical 
Observatory [104]. The major arm of the T was oriented 
east-west and had a length of 1300 m (96a at 22.25 MHz). The 
minor arm extended 312 m north from the center of the 
east-west arm. A photograph of the junction point, facing 
west down the east-west arm and showing the north-south arm 
going off to the right, appears in Figure 5.1. 

The array consisted of 624 full-wave dipoles strung 
between wooden poles. A reflecting ground screen of 
additional wires was situated a distance )/8 below the 
dipoles. 

The basic array element consisted of a pair of 
full-wave dipoles as pictured in Figure 5.2. For the north 


polar sky synthesis project [2], sets of four such basic 
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Figure 5.1. A View of the DRAO 
22.25 MHz Telescope 
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elements in the east-west arm were combined together to form 
48 sub-arrays, each 2d by 24, as shown in Figure 5.3. An 
additive combining network plus precise cable lengths for 
each of the four basic elements allowed control of the 
north-south pointing and beam shape for each sub-array. 

For the polar synthesis telescope the sub-array beams 
were centered on the north celestial pole. The 3 dB beam 
width in both the E-plane (east-west) and H-plane 
(north-south) was about 24 degrees, and the beam width 
between first nulls was about 60 degrees. 

For the interference-excising project two of the 
sub-arrays near the center of the east-west arm and a 
distance 350 m (26)). apart were chosen to form an 
interferometer. 

As the sub-array beam pattern had an approximate nul] 
for Cassiopeia A during transit, fringes from Cassiopeia A 
were initially not very strong. The antennas were changed 
part way through the observations. Instead of using the 
sub-arrays, signals were taken from just one of the four 
dipole-pair elements. A single element provided an antenna 
pattern centered on the zenith with an E-plane 3 dB beam 
width of 24 degrees and an H-plane width of 96 degrees. 


Signals from Cassiopeia A were then much improved. 
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FIGURE 5.3 Sub-arrays of Polar Synthesis Telescope 
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5.3 Cable Lengths and Losses 

One of the requirements for an interferometer is 
accurate Knowledge of the phases of signals from the 
antennas. To determine phase, the electrical lengths of all 
cables in a system must be Known. 

The only unknown cable lengths in the project (and also 
the longest) were for the cables pamneea thet ps tele 
amplifiers (located near the antennas) and the second mixers 
(located in the observatory building). Measurement of the 
cable lengths was assisted by the availability of a number 
of cable runs between the building and the antennas. Three 
cables, all of unknown length, were connected at the antenna 
end via a resistive Y splitter as in Figure 5.4. A signal 
generator transmitted an accurately Known frequency down one 
of the cables. Measurement of the amplitudes and relative 
phases of the signals at the terminated ends of all three 
cables were made. The generator frequency was then changed 
slightly and the measurements repeated. The signal generator 
was connected to the other two cables in turn and the above 
process repeated. 

From the amplitude measurements it was possible to 
determine the loss of each cable individually, and from the 
phase measurements at different frequencies the electrical 
lengths of each cable could be calculated. The length 
calculations were similar to those for a cable measurement 
procedure Known as the "extended probe" method [2]. 


The losses and electrical lengths at 5.0 MHz measured 
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for the cables to the east and west receivers are given in 
Table 5.1. The difference in losses is due to the use of 
low-loss heliax for a large part of the cable to the east 


receiver, rather than RG-8 as used for the west receiver. 


Table 5.1. Losses and Electrical Lengths of Receiver Cables 


at 5.0 MHz 
Receiver Cable Loss Electrical 
Length 
West 12.0 dB a} adel ON 
East 8.4 dB Jel BY 


5.4 Bandwidth Decorre lation 

A possible problem which may occur with radio 
telescopes is bandwidth decorrelation. If there is a large 
relative delay between two signals being correlated, a phase 
shift results across the bandwidth of the signals, thereby 
reducing sensitivity. 

Because the 52 KHz bandwidth of the FFT processor is 
small, bandwidth decorrelation is not serious. At the first 
null of the antennas in the east-west direction the position 
of a source would be 8 = 30 degrees from the zenith. With a 
baseline of D, = 26 wavelengths, the delay T) between the 


signal’s arrival at the two antennas is 
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D,. sin 9 


1 v 


=4 26 sin 30° as =f 
RiDIRIDRTase 0.59x10 sec 


An additional delay Ts due to the difference of 0.42 in the 
lengths of the IF cables should be added. 


: 0.4 
i Teo aera tiice i | 081 553) 1) rol ere 
5x108/sec 


The phase shift ¢ across a bandwidth of BW = 52 kHz for a 
delay of +t = 0.67x10-& seconds will be 
6° =*297- BW 


= 24x52x103/sec x 0.67x107© sec « 180/17 


ira 


This small phase shift will not produce any significant 
bandwidth decorrelation over the field of view of the 
antennas. 

If bandwidth decorrelation were a problem, it could be 
corrected for by adding a compensating amount of phase shift 
to the in-phase and quadrature spectra. The amount of 
compensation would be calculated to produce zero 


decorrelation for a desired point in the shy. 


5.5 System Operation 
The component parts of the entire system, including 
antennas, receivers, FFT processor and microcomputer were 


interconnected as outlined in the previous chapters. Figure 
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5.5 shows the equipment in the observatory building. The two 
second mixers and IF amplifiers were rack mounted above and 
below the Fluke synthesizer at right center of the photo. 
The FFT processor and microcomputer were enclosed in a 

_ shielding cage (end removed) to the left. A view of the FFT 
processor from above in Figure 5.6 shows the circuit boards 
in a rack at the back and the keyboard and readout console 
at the front. 

For the receivers the major operating considerations 
were frequency response and gain adjustment. As mentioned 
previously, the frequency responses of the second IF 
amplifiers were found to change due to mechanical 
instability of the IF transformer cores. Frequent monitoring 
and adjustment was therefore necessary. Receiver gains were 
set to provide a noise level of between 30 and 60 mVrms into 
the A/D converters. This level ensured that signal levels 
were above roundoff noise, even in the presence of strong 
sinusoidal interference. The largest sinusoid which would 
not cause significant clipping was about 300 mVrms. 

Most aspects of the system's operation were controlled 
by the microcomputer. For example, the microcomputer 
controlled the second local oscillator frequency and the 
integration time. For many of the observations the second LO 
was sequenced through four values: 4.483 MHz, 4.533 MHz, 
4.583 MHz and 4.633 MHz. The resulting receiver center 
frequencies were 22.325 MHz, 22.275 MHz, 22.225 MHz and 
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Figure 5.5. 
Equipment Layout in 


the Observatory Building 


Figure 5.6. 


The FFT Processor 








A major function of the microcomputer was the recording 
of observation results. In order for detailed analysis of 
the observations to be carried out at a later time, the raw 
spectra from the FFT processor were recorded on magnetic 
tape. Spectra directly from the FFT processor consisted of 
48 bit integers. These were converted by the microcomputer 
to 32 bit floating point numbers with a sign bit, a 7 bit 
twos complement exponent and a 24 bit unsigned mantissa. To 
reduce the amount of data stored on tape, spectra were 
recorded as 16 bit floating point numbers by truncating the 
mantissa to 8 bits. A header containing an identifying file 
number, the time of recording, the accumulation time, the 
local oscillator frequency and the number of clips counted 
during A/D conversion was recorded with each spectrum. A 
checksum for error detection was also included. 

The microcomputer employed a 40 rejection procedure for 
robust estimation of the centers of the in-phase and 
quadrature spectra. Estimation results were recorded on a 
two-pen chart recorder so that fringes could be observed as 
they were received. A second multiple-input chart recorder 
was used for the number of points deleted during robust 
estimation, the local oscillator frequency and the amount of 
clipping. An assembly language listing of the microcomputer 
program used during the observations is included in 
Appendix 3. 

A problem encountered with the microcomputer during the 


observations was eventually traced to the Am9511 Arithmetic 
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Processing Unit. This device was discovered on occasion to 
refuse to clear its BUSY flag. As a result, the 
microcomputer would get hung up in a loop waiting for 
arithmetic computations. A program modification which timed 
Am9511 operations and reissued commands if the device 
remained BUSY for too long was partially successful in 
overcoming this problem. However, further testing revealed 
that the device was operating incorrectly in other ways. 
Correct operation was found to be highly dependent on clock 
frequency which the manufacturer claimed could be anything 
from 0.3 to 3 MHz but which actually exhibited only a few 
windows a few tens of cycles wide near 3 MHz within which 
the device would work properly. The windows were dependent 
on temperature and would slowly drift, causing the author 
considerable consternation. 

This strange behavior appeared to be the result of a 
manufacturing defect. A second device (from the same 
manufacturing lot) performed in a similar manner. This 
integrated circuit was a relatively new and highly complex 
LSI chip, so manufacturing difficulties were not unexpected. 

The problem was finally alleviated by the simple 
addition of a heat sink to the IC. Though not called for in 
the manufacturer’s information, the heat sink Kept the IC's 
temperature sufficiently stable so that it would continue to 
operate satisfactorily once the clock frequency was set. 
Fortunately, spectra stored on tape were not affected by the 


difficulties with the arithmetic unit. 
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5.6 Data Analysis 

The spectra from the tapes were analyzed on the 
University of Alberta’s computing facilities. The FORTRAN 
program used for robust estimation is listed in Appendix 4. 
Analysis results presented in the next chapter were plotted 
by the computer on a CalComp plotter. 

Prior to robust estimation, a correction was made to 
the spectra to force the frequency response to be as flat as 
possible. This correction was necessary because of ripples 
in the frequency response of the receivers. 

An accurate reference for amplitude corrections was 
provided by the auto spectra. Auto spectra containing no 
interference were chosen, and fifth order polynomials fitted 
to estimate the receiver amplitude responses. All of the 
auto spectra received on the same night were then divided by 
these reference polynomials to produce a flat frequency 
response. 

The cross spectra were divided by the square root of 
the product of the polynomials for the two receivers in 
order to correct for amplitude ripples. Unfortunately a good 
reference for the relative phases of the two receivers was 
not available. Therefore, no correction could be made for 
differences in the phase responses of the receivers. 

The amplitude ripples found in the spectra before 
correction ranged from 0.15 dB to 2.4 dB, with the median 
amount of ripple being 0.88 dB. Spectra of strong fringes 


were examined for evidence of phase ripples, but no 
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detectable ripples were found. 
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6. Observation Results 
This chapter describes observation results obtained using 
the interference excising correlator at the Dominion Radio 
Astrophysical Observatory in Penticton, British Columbia 
during the period from November 12, 1979 to January 6, 1980. 
Section 6.1 contains a number of plots of in-phase and 
quadrature fringes both before and after interference has 
been removed. Some examples of raw spectra and their 
histograms are included. A number of plots of morning and 
evening observations are given, illustrating the range of 
interference levels encountered. 

Some of the fringes are plotted in terms of amplitude 
and phase in Section 6.2 in order to demonstrate the 
presence of scintillation. Variations in the amplitude of 
the fringes are compared to similar variations noted in the 
auto spectra. The performance of the first decile as an 
estimator for the auto spectra is investigated in Section 
6.5. The subsequent section describes an attempt at daytime 
observation employing an experimental frequency-changing 
scheme. 

An empirical probability distribution for interference 
amplitude is given in Section 6.5. The probability of 
interference is found to exhibit a power-law relationship to 
interference level. 

The chapter concludes with some general comments on the 


observations. 
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6.1 Plots of the Observations 

6.1.1 Examples with No Interference Removed 

Figures 6.1 and 6.2 show records of fringes from which 
interference has not been removed (i.e. the sample mean is 
used as the estimate of the centers of the co and quadrature 
spectra). In 6.1 fringes from Cassiopeia A (lower transit 
time 6:441, fringe per iod 17.0 min. at transit) are evident 
until 6:25 when interference obliterates them. Figure 6.2 
contains small fringes from Jupiter (transit time 4:00, 
fringe period 8.9 min.) with high levels of terrestrial 
interference at 5:30 and after 7:00. The peaks are truncated 
at +100, but some of those shown were actually as large as 
+1000. The interference peaks are periodic because the 
system is scanning sequentially through four adjacent 50 KHz 
regions of the spectrum between 22.15 and 22.35 MHz. Records 
of terrestrial interference such as these are commonly 


encountered in low frequency radio astronomy. 


‘All times in this chapter are Pacific Standard Time. 
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November 29, 1979 
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Figure 6.2. Fringes with Interference 


January 4, 1980 
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6.1.2 Examples with Interference Removed 

The same records as above are shown in Figures 6.3 and 
6.4 following interference removal via the robust estimation 
method of Section 4.6.3. Four traces are shown in these 
plots. Traces A and B are the in-phase and quadrature 
fringes. The solid lines are the robust estimates, while the 
x's are the sample means (no interference removed). The 
central regions of these plots are linear, whereas beyond 
+10 the plots become logarithmic in order to show the actual 
magnitude of the fringes when contaminated with 
interference. It should be remembered that the x’s represent 
the means of 108 spectral components, only a few of which 
contain interference. Therefore, in order to cause the mean 
to have a magnitude of 10° the interference, if it is 
contained predominantly in a single channel (as is usual), 
must have a magnitude of about 105 and is thus at least 
40 dB larger than the fringes (which are less than +10). 
Figure 6.4 contains some examples of interference of this 
amplitude being successfully rejected from the spectra. 

Trace C shows the expected standard deviation O,. of the 
under lying normal distribution as derived from the auto 
spectra using the estimation technique of Section 4.6.2. An 
approximate idea of the expected error in the fringes due to 
fluctuation noise can be obtained by dividing the value of 
trace C by Jin) (EIT), where n is the number of spectral 


components containing no interference and EIT gives the loss 


178 


, ot ae | a — 








Caer 
AY 
bevons oF conanetne 
























= = 


@ 2enugri af nworle ow ‘evods. as e a 


teudo7 seat etv tinsel ans sin ger wx 
ee ae 
aor) cit rwiode Sts 2a281 Wo? BOR at Jase 
oye Sat 
“ori 


s7u7 6 1bsup bra sas ] haaatial 
“Seudo1 ei a5 | gent tee 


‘es 


“rit 


orl 


am 


r 
a» © wre 
,co. Bitti ls 

* 


.(bsvemen sonet1Se?Asini on} ansem signee. 


foysd ese1sriw , seni. e186 atdia ezent 40 anol ge pee 
a: ¥ wie oF aeb10 ni iinet? t ager emosd ‘ate iq 6 3 

i?tw beisn sinoo. nerw zogrntay ait to: sbut 
| 213 tert? betgdmemes 4d bl worta rt eon 


viru Fo wet s yl ,BenenoqNas Isrt9eq2) aor to eae 
6! »xebwo af ,sretectenT panei 
msint ertt tor to obutt EIT en evs 


on pehes 
wey 2f 2 Lennano al on he & “ni vi tn tm ee he. 


’ 


279 


; ini aa 
fasai ts gurls vis @O7 feds. 19 aoudieg : 
; Z ’ Lau 17 a al a 
if: medd sas) ene fiotaws celia pets 


‘ty to sonstetaedars to 2st bail 


id =. 


.6%T and man 


19Q2 
tO .0 nhoOtTsives prebris ie 
' # "s 


sus 


ent. man} sav vo at é 


ea 
~ 


rs ue Le ; 

iS.¢ te notice? te" da at on votre 
- sas i i 4 

aeons }. oft a o ane t -< oA, 4 ® 


f Pe’ se ” 
’ iy 5 
5 ' P ‘ % 4 ' ; ¥ 
7 PREV Tg. 4 / : } — 
a = 


B) 


IN-PHASE FRINGES 


QUADRATURE FRINGES 


STO DEV 


DERETEO 


Ze 


Wake 


10 
tl 


de A AAARA LHL 
F ARS eeaes! 


Tos 
107 


LAMAR NT 
oy ee ee 





: . 5:30 6:15 7:00 
He aoe TIME 3 


Figure 6.3. Fringes after Interference Removal 


November 29, 1979 
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Figure 6.4. Fringes after Interference Removal 
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in stability of the estimate due to windowing. Trace D 
represents the percentage of components deleted by the 
robust estimation procedure (i.e. beyond 15.6. from the 
final estimate). This trace may be thought of as showing the 
percentage of points positively identified as outliers. 
Close-in contamination is not included though its effects 
are still minimized by the biweight estimator. 

From Table 4.2 the efficiency of the biweight estimator 
relative to the mean from a normal distribution is 1.041. 
The expected standard deviation oF of the fringes due to 


fluctuation noise will then be 


12 


; « aocmpeece (108) (EIT) ; EOres) 


For example, if 4Deleted=0 and EIT=0.44 then 0, =.148¢. wern 
most of Figure 6.3 only small percentages of the points are 
deleted, thus 0, =(0.148)(1.0)=0. 148 which is very smal] 
compared to the magnitude of the fringes. Similarly in 


Figure 6.4, 4, =0.53 which is small relative to the fringes. 
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6.1.3 Examples of Spectra and their Histograms 

Three examples of in-phase cross spectra and their 
histograms are given in Figures 6.5, 6.6 and 6.7. The 
examples correspond to cases of no interference, low-level 
interference, and strong interference. Each figure shows 
four consecutive in-phase spectra and their histograms. The 
spectra are all from observations on November 28-29 in 
Figure 6.3. The spectra and histograms are linear between 
the dotted lines but become logarithmic beyond the dotted 
lines in order to show the complete range of points 
necessary. 

Figure 6.5 shows four spectra from between 5:09 and 
5:13. In each case the histograms have a Gaussian shape with 
a standard deviation o, of about 1, which is as expected 
from the levels of the auto spectra (trace C in Figure 6.3). 
The centers of the four spectra change because the fringe 
from Cassiopeia A at this time is just passing through zero. 

Low to moderate interference appears in the spectra of 
Figure 6.6, taken from between 6:24 and 6:28. A central 
Gaussian distribution (o =1.3) is present with a few 
outlying points. 

Large amounts of interference are seen in Figure 6.7 
from between 7:15 and 7:19. The underlying Gaussian 
distribution eh SS), has been badly contaminated with 
interference. Between 25 and 45% of points are identified as 


outliers by the estimation prodedure for these spectra 
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(trace D in Figure 6.3), but there are sufficient numbers of 
non-contaminated points for the centers of the underlying 


normal distributions to be accurately determined. 
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Figure 6.5. Spectra and Histograms 


November 29, 1979 
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Figure 6.6. Spectra and Histograms 


November 29, 1979 
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Figure 6.7. Spectra and Histograms 


November 29, 1979 
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6.1.4 Morning Observations 

Figures 6.3 and 6.4 were both examples of morning 
observations, with increasing interference at the approach 
of sunrise. Figure 6.8 is another example, in this case 
continuing until after sunrise.. Fringes from Cassiopeia A 
(lower transit time 6:48) can be seen. A rapid increase in 
both levels and numbers of interfering signals begins at 
about 6:00, as would be expected due to the rapid increase 
in ionospheric electron density at sunrise. 

At 7:30, something unexpected happens. The fringe 
estimates begin to swing wildly, the standard deviation 
estimates jump erratically, and up to 100% of the 
observations are deleted. The estimates have broken down 
completely. The cause of this behavior was sweeping 
narrow-band signals which appeared over most of the 
shortwave spectrum at the same time nearly every morning 
during the observations. 

The sweeping signals were strong and easily observable 
on a spectrum analyzer. They continued throughout the day 
and gradually diminished as evening approached, usually 
becoming undetectable by about 16:00. Two or three separate 
sweeping signals were often seen simultaneously. Some 
started above 30 MHz and swept slowly down to below 15 MHz, 
while others appeared and disappeared rapidly and moved 
about the spectrum very unpredictably. The most troublesome 


swept signal at 22 MHz was one which recurred regularly with 
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Figure 6.8. Fringes after Interference Removal 


November 27-28, 1979 
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a repetition period of 10 seconds. 

Enquiries were made to the local Department of 
Communications office as to the source and legality of such 
sweeping signals. Monitoring stations eventually identified 
two locations producing the signals: Denver, Colorado and 
San Francisco, California. The signals were believed to be 
ionospheric soundings which are conducted worldwide for 
ionospheric research and radio propagation studies. A 
discussion of ionospheric sounding is contained in 
Davies [17]. Local Department of Communications officials 
indicated that the observed sweeping signals might be 
illegal, though confirmation of this has not been received. 

The presence of sweeping signals made observations 
beyond 7:30 and during the day impossible. Before 7:30, 
however, normal interference was excised successfully in all 
cases, as illustrated by the previous examples. Generally, 
the removal of interference consistently allowed from 60 to 
90 minutes (occasionally much more) additional observing 


time in the mornings during the field testing period. 
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6.1.5 Evening Observations 

In the evenings it was necessary to wait until some 
time after sunset before observations could begin. Strong 
interference caused the receivers to overload until 
ionospheric electron densities decayed sufficiently to 
reduce interference levels. After interference decreased, 
heavy scintillation was observed for a number of hours. 
Scintillation was apparent in all of the fringes observed 
during the field trials, and was always severe in the 


evenings. 


6.1.5.1 November 28 

One record of evening fringes from Cassiopeia A is 
shown in Figure 6.9. Transit time is 18:48 and the fringe 
period at transit is 17.0 minutes. A large amount of 
interference is initially present, with the strengths and 
numbers of interfering singals decreasing slowly over 1 1/2 
hours. Here again the received frequency is being 
sequentially scanned through four adjacent 50 KHz regions of 
the spectrum between 22.15 and 22.35 MHz. From the first 
half of trace D it may be noted that one of the four regions 
regularly has little or no interference present while the 
others have much more. The region of least interference is 
between 22.15 and 22.2 Mhz. 

From trace C the expected standard deviation varies 


from 19 to 6, implying a standard deviation of the fringes 
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due to fluctuation noise of from 2.9 to 0.9 or slightly more 
if a signicant number of points are deleted. Most of the 
scatter in the fringes is therefore due to fluctuation 
noise. The fringes are relatively noisy because the antenna 
pattern at this time is centered on the north celestial pole 
and has an approximate null for Cassiopeia A at upper 


transit. 
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Figure 6.9. Fringes after Interference Removal 


November 28, 1979 
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6.1.5.2 December 4 

A second evening record appears in Figure 6.10. The 
antenna pattern has been changed to have a maximum at the 
zenith, thus the signal from Cassiopeia A (transit time 
18:17) is much stronger as evidenced by an expected standard 
deviation of o =0.8, or fluctuation noise of o-=0.12. The 
receiving frequency is constant for this record and is 
centered on 22.225 MHz. 

A relatively constant amount of interference is 
received between 19:00 and 19:50, after which reception 
becomes quiet. There are many spikes in the fringes which 
are too large to be due to fluctuation noise. These appear 
to be due to strong scintillation and will be discussed 


further. in Section 6.2.1. 
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Figure 6.10. Fringes after Interference Removal 


December 4, 1979 


194 


d the a 
‘ a 25 
= Me 
iis t 
} 
he : 
tf 





it 








0 * “e 
othe OAS Qe Eh 
' 


ote 0s oaroe 


3 Le eee ar 
- eee 
> hr 


me oy are | 







6.1.5.3 December 5 

A record showing severe evening interference is shown 
in Figure 6.11. Very strong interference approaching 105 in 
some cases with up to 75% of points deleted continues from 
18:15 to 20:00. A single interferer would have to be 60 dB 
larger than the fringes to produce the levels shown. A 
considerable amount of clipping at the A/D converters was 
recorded a number of times until 19:45. 

The receiver frequency was again being scanned across 
four frequencies between 22:15 and 22:35 MHz. For most of 
the record, trace D shows one frequency which is usually 
free of any interference (22.15 to 22.2 MHz). 

The fringes from Cassiopeia A (transit time 18:13) are 
virtually unrecognizable for most of the record. For the 
first half of the record the fringes are being distorted 
because high interference levels are exceeding the dynamic 
range of the system. Periods of excessive clipping produce 
harmonics and contaminate the entire spectrum. Large single 
interferers such as those 60 dB above the fringe level 
produce spurious sidelobes (a result of FFT coefficient 
quantization) and thereby contaminate the rest of the 
spectrum. Finally, intermodulation distortion in the 
receivers was observed directly on a spectrum analyzer 
monitoring the first I.F. outputs during the time this 


record was made. Two very strong, steady shortwave signals 


were observed below 22.15 MHz which were clearly over loading 
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ene Tirst mixers and producing, intermodulation! products 
BeeorgjmiOmmeanc™ 1 irst lors pass*band. Wt swasson ly jatLem W9.30 
that the level of these two signals decreased 
sufficiently to allow normal observations to begin. 

After 20:00 the interference has become quite small. 
However, much larger variations are present in the fringes 
than would be expected due to fluctuation noise. These 


variations appear to be due to scintillation. 
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Figure 6.11. Fringes after Interference Removal 


December 5, 1979 
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6.1.5.4 December 6 

A final record of evening observations is shown in 
Figure 6.12. The received frequency was centered on 
22.175 MHz. Considerable interference occurs from 18:06 
until 19:15 with 10% to 30% (and in one instance 70%) of 
points being deleted. No clipping at the A/D converters was 
recorded during this record, but from 18:12 to 18:36 there 
is an increase in the levels of the auto spectra, as 
evidenced by trace C, which may be due to intermodulation or 
some form of broadband interference. | 

Fringes from Cassiopeia A (transit time 18:09) appear 
after 19:15. However, these fringes are badly distorted by 
scintillation. At a few points, such as at 19:12, large 
peaks occur. Examinations of the spectra for these points do 
not show any signs of the peaks being due to interference, 
hence they seem to be a result of severe scintillation and 
possibly momentary focussing of the radiation from 
Cassiopeia A on the antennas by the ionosphere. 

An interesting observation for this record is the 
extreme rapidity with which the interference disappeared. In 
a space of 20 minutes interference goes from occupying 40% 
of the spectrum to zero percent. Also, on this night the 
interference remained at zero from this time until 5:30 the 
nextamorning. scintillation continued until] at least 21:00. 
The absence of interference during the night was unusual and 


was seen only twice (December 5-6 and December 6-7) during 
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Figure 6.12. Fringes after Interference Removal 


December 6, 1979 
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GeemoCia tin ation 

6.2.1 December 4 

In Figure 6.13 the in-phase and quadrature fringes from 
part of Figure 6.10 are plotted in terms of amplitude and 
phase rather than real and imaginary components. The phase 
changes in a linear manner, as expected, but the amplitude 
fluctuates considerably. 

Nearly identical amplitude fluctuations can be observed 
riding on top of both auto spectra as illustrated in Figure 
6.14. Traces A and B are of the two auto spectra, showing 
only the variations on the tops of the spectra. In trace D 
the amplitude of the cross spectra from Figure 6.13 is 
redrawn. Trace C is most interesting, as it shows the 
geometric mean of the fluctuations in A and B. There is a 
striking similarity between C and D; in fact, they are 
almost identical. The same pattern of fluctuations is 
observed on the auto spectra. 

The conclusion is that a single broadband source with a 
rapidly varying amplitude is being received by both 
antennas. The source is Cassiopeia A and the amplitude 
variations are due to scintillation. According to the 
presently accepted view of scintillation as described by 
Briggs [23], irregularities in the ionosphere impose random 
phase variations on an incoming wavefront. At the earth’s 
surface, these phase variations are converted by 


constructive and destructive interference into a combination 
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of both phase and amplitude variations. Evidently, the 350 
meter baseline used for the above observations is short 
enough that phase variations appear nearly simultaneously at 
both antennas, as do the amplitude variations in Figure 
6.14. Therefore, the phase difference in Figure 6.13 changes 


in a fairly linear fashion. 
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Figure 6.13. Fringe Phase and Amplitude 
December 4, 1979 
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Figure 6.14. Auto Spectra, their’ Geometric Mean, 
and Cross Spectrum Amplitude 
December 4, 1979 
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6.2.2 December 6 

Two more plots of fringe amplitude and phase appear in 
Figures 6.15 and 6.16. These plots correspond to two 
portions of Figure 6.12. In 6.15, both amplitude and phase 
variations are evident. This is also true in 6.16, where 
there is good correspondence between peaks in both the 
amplitude and phase fluctuations. 

The auto spectra for 6.16 and the geometric mean of 
their variations are shown in Figure 6.17. As before, there 
is a close similarity among all four traces in this figure. 
However, on very close examination it may be seen that the 
peaks in trace B generally lead those in trace A by a very 
short time (approximately 20 seconds). The difference in the 
timing of the peaks is caused by ionospheric irregularities 
drifting between Cassiopeia A and the antennas with a 
velocity of at least 350m/20sec=17.5m/sec from west to east. 
The drifting of ionospheric irregularities has been observed 
in a similar manner since studies of the structure of the 


ionosphere first began (e.g. Booker [105], Briggs [23]). 
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Figure 6.15. Fringe Phase and Amplitude 
December 6, 1979 
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Figure 6.16. Fringe Phase and Amplitude 
December 6, 1979 
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Figure 6.17. Auto Spectra, their Geometric Mean, 
and Cross Spectrum Amplitude 


December 6, 1979 
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6.2.3 November 28 


A fourth record of fringe amplitude and phase, in this 
case corresponding to the fringes in Figure 6.3, appears in 
Figure 6.18. These fringes are from the lower transit of 
Cassiopeia A during the morning. Once again there are 
amplitude variations while the phase is quite linear. 
Scintillation is therefore occurring in the morning when 
ionospheric electron density should be at a minimum. As 
Cassiopeia A is just above the horizon at this time 
scintillation is increased by the long path the signals must 


take through the ionsphere. 
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Figure 6.18. Fringe Phase and Amplitude 
November 29, 1979 
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6.3 Estimation for the Auto Spectra 

Two examples demonstrating the performance of the first 
decile as an estimator for the location parameter of an auto 
spectrum are presented in this section. 

Figure 6.19 shows estimates for one of the auto spectra 
associated with the fringes in Figure 6.10. In trace A, the 
solid line represents an estimate derived from a +3, 
rejection procedure. The x’s are uncorrected estimates of 
) 


Paty: 
Trace B gives the standard deviation of the points remaining 


the first decile (i.e. the eleventh order statistic, 


after iterative +30 rejection, and C gives the percentage of 
points deleted. 

A comparison of the solid line and the x’s in trace A 
indicates that the first decile is consistently below the 
rejection estimate by about 2.0, or 2.0/100=2%. The 
integration time during this record is 20 seconds, which 
corresponds to K=8137. The correction factor required in 
equation 4.63 is then (1+1.3/¥8137)=1.014. After correction, 
the first decile estimates of O, will therefore be 
(1-0.02)(1.014)=0.994 times the 30 rejection estimates, 
producing a bias of -0.6%. The slight discrepancy is due to 
a small amount of ripple in the auto spectrum. 

The expected standard deviation of the auto spectrum 
from equation 4.58 is 100//8137=1.11. From trace B, the 
standard deviation after deletions is actually about 1.5, as 


would result if a ripple of something less than 1% were 
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Figure 6.19. Auto Spectrum Estimates 


December 6, 1979 
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present in the auto spectrum. The fluctuations in the level 
of the auto spectrum are due to scintillation of Cassiopeia 
A, as shown in Figure 6.13. 

A second example of auto spectrum estimation is given 
in Figure 6.20, corresponding to the fringes in 6.4. In this 
example, the first decile estimates are less than the 
rejection estimates by about 10.0, or 10.0/830=1.2%. The 
correction factor with K=24411 is (1+1.3//24411)=1.0083. The 
final estimate is then (1-0.012)(1.0083)=0.996 times the 
correct value, resulting in an error of -0.4%. 

Overall, the first decile performs very well as an 


estimator for the centers of the auto spectra. 
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Figure 6.20. Auto Spectrum Estimates 
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6.4 An Attempt at Daytime Observation 

From an interference standpoint, daytime in the winter 
of a sunspot maximum is the worst possible time to attempt 
low frequency astronomy. Indeed, the levels of interference 
were found to be very high in the daytime. Intermodulation 
distortion in the receivers initially made daytime 
observations impossible. Modifications to the RF and first 
mixer stages of the receivers detailed in Chapter 3 were 
undertaken to improve their intermodulation performance. 

An attempt at daytime observation was made on 
January 6, 1980. The major objective of this attempt was to 
test a scheme of automatically changing the receiver center 
frequency if excessive interference was encountered in order 
to look for a quieter part of the spectrum. 

A fortunate circumstance was the absence of sweeping 
interference on the day of the test, probably because the 
day was a Sunday. Sweeping interference was not detectable 
on a spectrum analyzer at any time during the day, as it had 
invariably been during the days of previous observations 


which were all weekdays. 


6.4.1 Frequency Changing 

A very simple method of frequency changing was 
employed. The spectrum between 22:15 and 22.35 MHz was 
divided into four 50 KHz slots, any of which could be 


selected for observing. Changes were always made in order, 
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starting from the highest slot, 22.325 MHz, through 22.275, 
22.225 and 22.175 MHz and then back to 22.325 MHz. 

The center frequency was changed on either of two 
conditions; 1) an excessive amount of clipping occurred at 
the A/D converters due to very large interference; 2) more 
than 40% of spectral components were deleted due to a large 
number of small interferers. The above scheme was easy to 
implement and seemed adequate to test the merits of having a 


frequency-changing capability. 


6.4.2 Daytime Results 

The results of the daytime test were mixed. 
Interference levels were high enough that, even with the 
improved receivers, serious intermodulation was at times 
noticeable on the spectrum analyzer. In addition, 
interference of a broadband nature (100 kHz or more) was 
sometimes seen. However, there were also times when two of 
the 50 KHz slots, centered on 22.175 and 29.295 MHz, were 
virtually free of interference. 

Some fringes from Cygnus A (transit time 12:56) were 
seen between 10:00 and 14:00, but these were interspersed 
with periods of severe interference which the system could 
not remove completely. The fringes were extremely irregular 
in both amplitude and phase. A portion of the daytime 
observations can be seen in Figure 6.21. 

The ability to change frequencies was of very limited 


Vollcueihestwounigheresiots, 22.275 and 22.325 3Mnz, 
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contained strong interference at all times. The system 
therefore spent more than 99% of its time in the two lower 
slots. The system changed between these two frequencies only 
three times during an eight hour observing period. 

Frequency changing would be of the greatest advantage 
if there were many regions of the spectrum from which one 
could choose, with an equal and independent probability of 
interference at any frequency. Such was not the case during 
the daytime test above, hence frequency changing was of 


little benefit. 
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6.5 A Probability Distribution for Interference 

A very important point in predicting the success of 
interference excising is some Knowledge of how interference 
tends to be distributed in amplitude. In particular it is 
desirable to Know to what degree very small and undectable 
interference can be expected to be present. 

An analysis of the records of the observations made at 
the Dominion Radio Astrophysical Observatory was conducted 
to determine whether or not the distribution of interference 
followed a pattern. The analysis consisted of producing 
histograms of the magnitude of deviations of cross spectral 
components from the centers of the spectra. The fringes 
found by robust estimation were subtracted from the co and 
quadrature spectra to center the spectra about zero, and 
then the rms deviation (Vco2+quadrature? ) of each component 
of the spectrum was found. The rms deviations of many hours 
of observing were tabulated to form each histogram. 

Because interference varied over a range of about 
60 dB, it was impractical to perform tabulation in a linear 
manner. Instead, a tabulation with equal sized intervals was 
Pe onned on the logarithms of the rms deviations. In 
effect, the size of the intervals or bins into which 
spectral components were collected then increased in 
proportion to the amplitude of the components. The number of 
counts in each bin was then divided by the size of the bin 


to restore linearity. 
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A histogram of results from part of the night of 
December 6-7 when almost no interference was received 
appears in Figure 6.22. The dashed line represents a 
Rayleigh distribution, which is the expected distribution 
with no interference. Only 0.031% of the 38,700 points 
tabulated in this histogram were found to deviate 
significantly from the Rayleigh distribution. 

A second histogram, this time with considerable 
interference, is shown in Figure 6.23. This histogram 
corresponds to the first two-thirds of the record in Figure 
6.10. There is now interference extending for over four 
decades past the end of the Rayleigh distribution. The most 
important feature is the nearly linear manner in which the 
histogram decreases over the four decades of interference. A 
power-law relationship between the probability of 
interference and its magnitude is suggested. 

Five additional histograms containing interference are 
presented. Figure 6.24 shows a period of considerable 
interference from the evening of January 4, while 6.25 is a 
period of very small interference later the same night. 
Figure 6.26 is from the daytime observations of January 6. 
Figures 6.27 and 6.28 are from the evening observations of 
December 5, but for different center frequencies (22.175 and 
22.325 MHz, respectively). For 6.28, spectra from periods of 
excessive clipping have been excluded to prevent distortion 
products from contaminating the histogram. 


In each case above, the probability of interference is 
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functionally related to amplitude by a power law. The slopes 
of the linear portions of the histograms vary from -1.23 to 
-1.33, with the average being -1.3. It therefore appears 
that the probability density function for interference has 


the form 
f£(p) = apts (Ooms) 


over the range of amplitudes shown in the histograms, where 
p is the interference power anda is a constant of 
proportionality which varies with the degree of interference 
being received. The above power law does not hold true 
beyond a normalized count of about 10-4 in most of the 
histograms, perhaps due to there being too few points to 
count accurately and also due to limits on the system’ s 
dynamic range. Also, the power law cannot hold as p 
approaches zero, as f(p) would go to infinity. 
Unfortunately, the behavior of f(p) as p becomes small 
is masked by the Rayleigh distribution of fluctuation noise. 
Because fluctuation noise is contained in all spectral 
components, the combined distribution is actually a 
convolution of the Rayleigh distribution and the 
interference distribution. The interference distribution 
alone, without fluctuation noise, includes a delta function 
at its center which represents components with zero 
interference. The very regular nature of the interference 
distribution for a number of decades beyond the Rayleigh 


distribution suggests that the power law may hold for 
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interference at lower levels as well. It is only 
contaminated points near the junction of the interference 
distribution and the Rayleigh distribution which are 
important, because smaller interference will not affect 
location estimates significantly and larger interference is 
rejected by the robust estimation procedure. The power law 
should perform reasonably well for predicting average. 
numbers of interferers near the junction point. 

It must be emphasized that the power law found above is 
derived from the long term behavior of interference (over a 
few hours, at least) and as such is only a statistical 
average. Over the short term, interference is less 
predictable and the power law loses much of its 
Significance. Also, the proportionality factor a in equation 
6.2 is not a constant, but will exhibit its own statistical 
distribution reflecting ionospheric activity. A 
determination of the dependence of a on such factors as the 
time of day, the season and the solar cycle would require a 
long term program of observations and data collection. 

A relative idea of the range of a encountered during 
the observations can be found by measuring the distance 
between the peak of the Rayleigh distribution and the 
junction point of the Rayleigh and linear sections of the 
curves in the histograms. The junction point varies from a 
factor of 10*° below the peak in Figure 6.24 to 10278 below 
inepigune 6.25, implying a change ino by a factor of 


10%8=63 between these records. Over this range of a the 
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Figure 6.22, Histogram of Cross Spectrum 
Probability vs. Amplitude 
December 6-7, 1979 
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Figure 6.23. Histogram of Cross Spectrum 
Probability vs. Amplitude 
December 4, 1979 
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Figure 6.24. Histogram of Cross Spectrum 
Probability vs. Amplitude 
January 4, 1980 
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Figure 6.25. Histogram of Cross Spectrum 


Probability vs. Amplitude 
January 47.1980 
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Figure 6.26. Histogram of Cross Spectrum 
Probability vs. Amplitude 
January 6. "L900 
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Figure 6.27, Histogram of Cross Spectrum 
Probability vs. Amplitude 
December 5, 1979 
Frequency, = 22.175 Miz 
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Figure 6.28. Histogram of Cross Spectrum 
Probability vs. Amplitude 
December 5, 1979 
Frequency = 22.325 MHz 
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power law f(p) =ap~*%3 remains unchanged. 

It must be noted that the value of a, or alternatively 
the incidence of interference, is strongly dependent upon 
the observing frequency. This fact is clearly demonstrated 
by a comparison of Figures 6.27 and 6.28, which were 
tabulated at the same times but at frequencies 150 kHz 
apart. The values of o are different by a factor of 


approximately 10 for these two figures. 


6.6 Comparison to the Log-Normal Distribution 

A previous study by Wheeler [31] found that the 
statistical distribution of interfering signal powers could 
be approximately described by a log-normal function. The 
interference distributions in the preceeding section also 
fit portions of log-normal curves reasonably well. The 
author stresses that the fitting of his own observations or 
others to a curve such as a log-normal or power law is 
purely empirical and should not be extrapolated or 
generalized to any degree without further evidence. 

The fitting of interference distributions to a 
log-normal function has one serious drawback, namely a 
comewiat arbitrary choice of function parameters. To 
illustrate this difficulty, two log-normal curves with 
different parameters are shown in Figure 6.29. The 
expression for a log-normal density function is given in 
equation 1.7. Over a range of p from 10-2 to 10° (ten 


decades) the two curves in Figure 6.29 are nearly identical, 
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with a slope (on the log-log plot) ranging from -1.0 to 
-2.0. A portion of either curve over a selected four decade 
range within this region would undoubtedly fit the observed 
data, which has a fairly linear slope of about -1.3. 
However, there are no unique choices of parameters m and o 
which produce the desired fit. 

Note that the two curves in Figure 6.29 diverge sharply 
as p becomes small. Measurements of interference levels down 
to extremely low values of p would be necessary to determine 
whether any particular distribution is truly representative. 
The author believes that the equipment designed for this 
interference excising project may be the most sensitive ever 
used for measurements of interfering signal levels. 
Wheeler’s method, for example, was not as sensitive because 


averaging to reduce background noise was not possible. 


6.7 General Comments on the Observations 

A number of general comments on the observations may be 
made. First of all, terrestrial interference during the 
observations was extremely common. In addition to mornings 
and evenings, a great deal of interference was received 
during the nights. The incidence of nighttime interference 
changed during the observing period. Between November 20 and 
November 30, low-level interference was virtually continuous 
each night. On December 5 and 6, the nights were free from 
interference. The remainder of the records showed intervals 


of low-level interference at night lasting from a few 
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minutes to many hours. 

The part of the spectrum chosen by the Dominion Radio 
Astrophysical Observatory for their 22.25 MHz telescope is 
extremely good from an interference standpoint. A view of 
daytime radio transmissions on a spectrum analyzer showed 
signals to be less common and generally lower in amplitude 
by as much as 30 dB or more near 22.25 MHz than in 
neighboring parts of the spectrum. Even within this region 
of low interference, some frequencies were better than 
others. As indicated by the results in this chapter, there is 
a marked difference in the incidence of interference across 
the 200 KHz band from 22.15 to 22.35 MHz. There was 
consistently less interference in the lower half of this 
band, and the 50 KHz region from 22.15 to 22.20 MHz 
generally had the least of all. 

Scintillation was continually present during the 
observations. As both the incidence of terrestrial 
interference and scintillation are strongly correlated with 
solar activity, the amounts of both experienced during the 
field testing were not unexpected. 

Overall, the interference excising system performed 
very well as long as interference levels were not too large. 
Generally this was the case for periods of from 60 to 90 
minutes in both the mornings and evenings when interference 
was increasing or decreasing. In addition numerous occasions 
of low-level interference were encountered during the night. 


Interference was removed successfully on these occasions. 
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7. Conclusions 
What are the conclusions to be reached following this 
project? First of all, interference detection and removal 
through spectral analysis and robust estimation has a great 
deal of promise for enhancing the quality of low frequency 
observations and for extending the amount of time during 
which observations may be conducted. Second, the very 
important aspect of the distribution of interference 
amplitudes requires more investigation before a definitive 
answer as to how much improvement is possible can be given. 
This chapter begins by outlining the performance of the 
present system and noting additions or changes which could 
be made. The question of the distribution of interference is 
then examined, and recommendations for further research are 


made. 


7.1 System Performance 

A number of comments about the receivers, the FFT 
processor and the robust estimation method may be made. 

Overall, the receivers worked well but could certainly 
be improved upon. In particular, dynamic range and frequency 
response stability were two problem areas. During the 
reception of strong interference, the avoidance of nonlinear 
operation in any stage of the receivers is crucial. A 
calibration system which would have allowed the monitoring 


of both amplitude and phase response of the system would 
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have been desirable for maintaining proper filter tuning and 
for allowing correction of residual response ripples. 

The FFT processor, though not designed for large 
interfering signals, was certainly adequate for the majority 
of signals encountered during the observations. Improvements 
in dynamic range could be made through the use of more than 
8 bits for FFT processing, but the benefits such as the 
incremental gain in observing time and quality should be 
weighed against the cost. An increase in spectral 
resolution, as discussed further on, would undoubtedly be 
beneficial. The calculation of auto spectra as well as cross 
spectra is quite valuable for robust estimation. 

Careful containment of electromagnetic interference 
generated by the digital circuitry is most essential. An 
oscilloscope probe connected to the FFT processor and 
brought outside the shielded enclosure was enough to radiate 
detectable RF interference. Electromagnetic interference 
produced by digital equipment at low frequencies could be a 
serious problem to radio astronomy. 

The robust estimation procedure developed for the 
project worked very well. Some possible improvements, such 
as the use of adaptive estimation, have already been 
mentioned. Additional possibilities could include the 
identification of specific modulation types and the deletion 
of an appropriate amount of bandwidth, or the use of past 
spectral information to determine if some frequencies should 


be weighted less heavily than others due to a higher 
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probability of interference. 

A major difficulty encountered during the observations 
was the occasional presence of signals which were not fixed 
in frequency or were not narrow in bandwidth. Sweeping 
signals from ionospheric sounding were one such problem. If 
sweeping signals could be readily identified ante WrRe 
present only intermittently, as they appeared to be during 
the field tests, then it should be possible to operate 
between the sweeps. A more serious potential problem could 
be spread-spectrum communications systems used by the 
military. The author does not Know to what extent such 
systems are presently in use or at what frequencies, but 
attempts should be made to find out. 

The interference power distribution observed during the 
field trials is quite interesting due both to its 
consistency (a power law f(p) =ap~*3) and its variability 
(a wide range of a depending upon frequency and time of day, 
and varying considerably from one day to the next). The 
relatively weak dependence of the probability of 
interference upon interference power (as compared to an 
exponential dependence, for example) suggests that if a is 
large enough for significant interference to be present, 
then the largest of the interfering signals should certainly 
be detectable via spectral analysis. Also, if the number of 
detectable interfering signals is not too large, then it is 
unlikely that there will be many just below the detection 


limit which could produce significant errors in observation 
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results. The obvious conclusion is that an interference 
detection and removal system as described in this thesis 
will be extremely successful if the number of interfering 
signals is small. If the amount of interference is too large 
then errors due to undetected interference will at some 
point begin to dominate. 

The results in this thesis show the interference 
excising system operating on occasion with up to 50% of 
channels being deleted because of interference. More 
typically, successful operation is possible with 20 to 30% 
of the overall bandwidth being rejected. The nominal 
increase in observing time without disruption from 
interference in the mornings and evenings is from 60 to 90 
minutes. In addition, night-time interference was noted on 


the majority of nights and was removed with no difficulty. 


7.2 Recommendations for Further Research 

This thesis has demonstrated that interference excising 
is highly successful during the winter nights of a solar 
maximum. However, it is not successful during the days. The 
main question which arises is how well such a scheme can 
operate during a solar minimum when conditions are more 
favorable for astronomy. As typical critical frequencies 
during the summer day of a solar minimum are nearly the same 
as those during the night of a solar maximum, it is 
conceivable that interference excising could allow very long 


periods of continuous, interference-free observing during 
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the summer and possibly the winter of a solar minimum. 
Further research should therefore include the use of the 
present or a similar system in conjunction with a large 
decametric telescope to evaluate its performance during the 
next solar minimum in 1984-86. 

As the probability of interference appears to be highly 
dependent on frequency, studies could be undertaken to 
identify other regions of the spectrum below 22.25 mHz where 
an interference excising system of the type described herein 
would be most effective. In particular, spectral regions 
with very narrowband signals and phase-coherent carriers 
should be sought. 

Far more about the probability distribution of 
interfering signals needs to be Known. Is it similar at all 
frequencies and for all classes of signals? How does it vary 
from day to day, season to season, and over the sunspot 
cycle? And, most importantly from an interference excising 
viewpoint, how is interference distributed below the 
threshold of detectability where it can produce observation 
errors? In order to answer this question, more sensitive 
measurements of interfering signals are required. The 
simplest method would be higher resolution spectral 
analysis, perhaps down to a few Hertz. The basic question to 
be asked is does the interference distribution turn over at 
some point (like the log-normal curve) and if so, at what 
point? 


There is a massive amount of literature on ionospheric 
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radio propagation which would shed a great deal of light 
upon and perhaps even answer some of the above questions. 
The decreasing costs of digital electronics and the 
improving performance of signal processing IC’s makes 
interference detection and removal quite feasible for 
decametric telescopes. The ability to avoid terrestrial 
interference using the methods developed in this thesis wil] 
undoubtedly be of much benefit to low frequency radio 


astronomy. 
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Appendix 1 - A Derivation of Equivalent Integration Time 


Consider two real Gaussian series f and g with zero means 
and variances of? and 09>, respectively. Assume that 


consecutive samples are independent, hence 


ov[ i f ee E Pete | ie >] af j 
c is ;) 1 J 
0 9 17 j 


= 2 i=) 
Ol ar 


let the covariance of f; and 9; be nonzero only if i=j, thus 


Cov[f,, g,] = PT.S» i=j 
0 » ij 


where pis the correlation coefficient of f and g. 


Define the product P of f and g to be 


Finding P is equivalent to correlating the two series. The 
expected value and variance of P can be found by considering 
P to be a linear function of a random variable fg. It is 
well Known that the expected value and variance of a 


weighted sum U of independent terms Y; are given as follows: 
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The relative error RE in P is given by the ratio of vV[P] 


and E[P]: RE, = /TIPT 





This result is as expected for a rectangular window where 
all points are weighted equally. 
Now consider two weighted series f’ and g’ formed by 


multiplying f and g by a set of window coefficients w;: 
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Equivalent integration time (EIT) for a windowed series 
may be defined as the square of the ratio of the relative 


error with no window to that with the window: 
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Some examples of equivalent integration time for a few 
windows using the above expression are EI7T=1.0 for a 
rectangular window, EIT=0.516 for a Hanning window, and 


EIT=0.444 for a Kaiser-Bessel window with parameter a=2.5. 
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Appendix 2 - Products of Gaussian Variables 
1. Expected Value of the Product of Two Gaussian Variables 


Start with two independent standard normal variables x, 
and X5,. Two partially correlated Gaussian variables y, and 
Ya may be created as follows: 

Veena lor 


Voi Bo, (pate ylepexyR O<p<l 


Then 


Viy,] = 0,7E[x,] = ane 
Ely,] = o,pE[x, ] + O, Lapce [xi] = 0 
My oles 5, *p-V[x, ] ~ 6p 21-p*)V [3x ] 
z 0,2 (p2+1-p?) i He 
The expected value of the product of y, and Yo is 
Elyyy,] = E[(0,x,) 95 (px, + v1-p*x,)] 
= Elo,0,px,° + Say aE eee 
0 
2 2 Lease 
= 05 05( Elx,*] + Y1-p*ELse%5] 


= Diners 
0, 5pE[x, ] 0155p 


as X, and Xp are independent. 
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2. Variance of the Product of Two Gaussian Variables 


The variance of the product of y, and Yo defined above 


will be given by 


= F- 2? 
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To find E[x,4], define the fourth moment as 
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Using integration by parts, let 


3 du = 3x_2 
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Then, using standard integration tables, 
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Finally, 


3. Expected Value of 


Start with four 


variables X41 Xo, XQ 
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the Product of Four Gaussian Variables 


independent standard normal N(0, 1) 


and Xa. It is desired to create four 


Gaussian variables F,, Fs, G, and Gj with the following 


variance-covariance matrix: 
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These variables will represent real and imaginary components 


of the complex spectra of two Gaussian series. The desired 


relationships above can be obtained by letting 
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With four variables as defined above, it is possible to 


find the expected value of the product of all four: 
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Appendix 3 - Microcomputer Observing Program 


The following is an MC6800 assembly language listing of the 
program used to control the microcomputer and FFT processor 
during observations. The program enables DMA transfers after 
a specified integration time, displays and records the 
spectra, and performs robust estimation using 4 standard 
deviation deletions. Estimation results, numbers of points 
deleted, amounts of clipping, and second local oscillator 


frequencies are recorded on chart recorders. 
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0000 1 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
000093 
00010 
0001 1 
00012 
00013 
00014 
00015 
00016 
00017 
00018 
00019 
00020 
0002 1 
00022 
00023 
00024 
00025 
00026 
00027 
00028 
00029 
00030 
0003 1 
00032 
00033 
00034 
00035 
00036 
00037 
00038 
00039 
00040 
0004 1 
00042 
00043 
00044 
00045 
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0040 
0040 
0042 
0044 
0046 
0048 
OO4A 
004B 
004D 
OO4F 
0050 


0084 
0088 
OO08C 
0o08D 
009 1 
0092 
0094 
0098 
OO39C 
OOAO 


0002 
0002 
0002 
0002 
0002 
000 1 
0002 
0002 
0001 
0001 
0001 
0002 
0001 
000 1 
000 1 
0001 


0004 
0004 
000 1 
0004 
000 1 
0002 
0004 


0004 
0002 


SCALE 
OFFSET 


NAM 
ORG 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 
RMB 


DISP 
$0040 


NAPHANHDH HHA HK HHBR HRB KBAAWWWHA WANA BAD ABAA A DNHANNNNAD 


Release 1.1 


TEMP INDEX STORAGE 
RESULTS POINTER 

DATA POINTER 

DELETION POINTER 

DISPLAY POINTER 

WORKING REGISTERS FOR 
FLOATING POINT CONVERSION 


CHANNEL POINTER 
INTEGRATION COUNTER 


FIRST FFT INDICATOR 

DELETION VECTOR POINTER 

COUNTER FOR NO. OF CHANNELS NOT DELETED 
FILE HEADER 

FILE NUMBER 

TIME OF DAY 

INTEGRATION COUNT 

SECOND L.O. FREQUENCY 

CLIPPING COUNTER 

INTEGRATION COUNTER 

CO SPECTRUM LOCATION ESTIMATE 

CO SPECTRUM NO. OF CHANNELS DELETED 
CO SPECTRUM STD. DEV. ESTIMATE 
QUADRATURE SPECTRUM PARAMETERS 


AUTO1 SPECTRUM PARAMETERS 


AUTO2 SPECTRUM PARAMETERS 


TAPE CHECKSUM 

TAPE STORAGE POINTER 
DELETION MAXIMUM 
MAX=FACTOR * STD. DEV. 
CRT DISPLAY PARAMETERS 
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DISP 


00046 
00047 
00048 
00049 
00050 
0005 1 
00052 
00053 
00054 
00055 
00056 
00057 
00058 
00059 
OOO060 
0006 1 
00062 
00063 
00064 
00065 
OOO66 
00067 
00068 
00069 
00070 
0007 1 
00072 
00073 
00074 
00075 
00076 
00077 
00078 
00079 
00080 
0008 1 
00082 
00083 
00084 
00085 
00086 
00087 
00088 
00089 
OO0O0sO0 
0009 1 
00092 
000393 
00094 
000395 


OOA2 
OOA3 
OOA4 
OOAG6 
OOA7 
OOA8 
OOAS 
OOAA 
OOAD 
OOBO 
OOB3 
OOB6 


0100 
0100 


0001 
0001 
0002 
0001 
0001 
0001 
0001 
0003 
0003 
0003 


000 1 


0080 


3227 
3370 
37AC 
3600 
3FOO 
3F80 
4000 


BAOO 
D202 
D204 
D206 
D208 
D20A 
D20C 
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PHASE RMB 1 ADD/SUBTRACT/SWITCH SPECIFIER 
EXP RMB 1 EXPONENT 
DPTR RMB 2 SPECTRUM POINTER 
RAMPH RMB 1 DISPLAY RAMP 
OLDN RMB 1 OLD NO. OF CHANNELS DELETED 
ITMAX RMB { MAX. NO. OF ITERATIONS 
TMRDEL RMB 1 ITERATION DELAY COUNTER 
NFREQ RMB iS} NEXT L.O. FREQUENCY 
FMAX RMB 3 MAX. L.O. FREQUENCY 
FMIN RMB 6) MIN. L.O. FREQUENCY 
FINC RMB 3 FREQUENCY INCREMENT 
DMAFLG RMB 1 DMA FLAG 
* 
ORG $0100 
WCOEF RMB 128 WINDOW COEFFICIENTS 
* 
FPSUB EQU $3227 FLOATING POINT SUBTRACION ROUTINE 
CLRMEM EQU $3370 SUBROUTINE TO CLEAR MEMORY LOCATIONS 
CRLF EQU $37AC CARRIAGE RETURN, LINEFEED ROUTINE 
FPOUT EQU $3600 FLOATING POINT TTY OUTPUT ROUTINE 
DVCTR1 EQU $3FOO DELETION VECTOR 1 
DVCTR2 EQU $3F80 DELETION VECTOR 2 
LATCH EQU $4000 OUTPUT LATCHES (NOT USED) 
TOS EQU $5000 TOP OF AM9511 STACK 
CMSTAT EQU $5001 AM9511 COMMAND AND STATUS REGISTER 
TIMER EQU $5008 MC6840.TIMER 
TTYPIA EQU $8004 TTY OUTPUT REGISTER 
RDRCTL EQU $8007 READER CONTROL REGISTER 
KEYBD EQU $8008 KEYBOARD REGISTER 
ACIACR EQU $8010 ACIA CONTROL REGISTER 
ACIATR EQU $8011 ACIA DATA REGISTER 
FCNTRL EQU $8020 FFT CONTROL PIA 
FFTFLG EQU $8021 CONTROL REGISTER WITH FFT DONE FLAG 
PCNTRL EQU $8022 FFT PHASE CONTROL PIA 
CLKFLG EQU $8023 CONTROL REGISTER WITH FFT CLOCK FLAG 
HZ EQU $8042 1 HZ CLOCK PIA 
WINPIA EQU $8600 WINDOW CONTROL PIA 
XTMP EQU $AOOF TEMPORARY INDEX STORGAE 
DCO EQU $BOOO0 CO SPECTRUM DISPLAY VALUES 
DQUAD EQU $B 100 QUADRATURE DISPLAY VALUES 
DAUTO1 EQU $B200 AUTO1 DISPLAY VALUES 
DAUTO2 EQU $B300 AUTO2 DISPLAY VALUES 
co EQU $B800 CO SPECTRUM VALUES 
QUAD EQU $BAOO QUADRATURE SPECTRUM VALUES 
DAC2 EQU $D202 7 D/A CONVERTER REGISTERS 
DAC3 EQU $D204 
DAC4 EQU $D206 
DAC5 EQU $D208 
DAC6& EQU $D20A 
DAC7 EQU $D20C 
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DISP 


OO0096 
00097 
00098 
O0O00gs9 
00100 
00101 
00102 
00103 
00104 
00105 
00106 
00107 
00108 
00109 
00110 
00111 
00112 
00113 
00114 
00115 
00116 
00117 
00118 
00119 
00120 
00121 
00122 
00123 
00124 
00125 
00126 
00127 
00128 
00129 
00130 
00131 
00132 
00133 
00134 
00135 
00136 
00137 
00138 
00139 
00140 
00141 
00142 
00143 
00144 
00145 


3000 


3000 
3003 
3006 
3008 
300A 
300C 
SOOE 
3010 
3012 
3014 
3016 
3018 
3SO1A 
301C 
SO1E 
3020 
3023 
3025 
3028 
302A 
302D 
302F 
3032 
3035 
3038 
3039 
303C 
303E 
3041 
3044 
3046 
3049 
304B 


3O04E 
304F 
3052 
3054 
3056 
3058 
305A 
305C 


Zoo 
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D20E DAC8 EQU $D20E 
E1D1 OUT EQU $E1D1 TTY CHARACTER OUTPUT ROUTINE 
* 
ORG $3000 
* LOAD WINDOW COEFFICIENTS FROM RAM 
* TO WINDOW COEFFICIENT MEMORY 
BD 3840 BEGIN USR INFREQ 
CE 8600 WINDOW LDX #WINPIA INITIALIZE WINDOW PIA 
86 04 LDA A #4 
C6 FF LDA B- #$FF 
6F O1 CLR 1,X 
6F O3 CLR 3,X 
E7 OO STA B X 
A7 O1 STA A 1,X 
86 FO LDA A #$FO 
A7 OO STA A X 
86 F2 LDA A #GF2 
A7 OO STA A X 
E7 02 STA B 2,X 
86 04 LDA A #4 
A7 O3 STA A 3,X 
CE OOFF LDX #WCOEF-1 
A6 ‘OO OVER LDA A xX 
B7 8602 STA A WINPIA+2 
86 F6 LDA A #$FE6 
B7 8600 STA A WINPIA 
86 F2 LDA A #$F2 
B7 8600 STA A WINPIA 
7C 8600 INC WINPIA 
7A 8600 DEC WINPIA 
08 INX 
8C 0180 CPX A#AWCOEF+128 
26 E5 BNE OVER 
7F 8603 CLR WINPIA+3 
7F 8602 CLR WINPIA+2 
86 FE LDA A #$FE 
B7 8600 STA A WINPIA 
86 04 LDA A #4 
B7 8603 STA A’ WINPIA+3 
* 
* INITIALIZE PERIPHERALS AND VARIABLES 
* 
OF INIT SEI 
CE 8004 LDX #TTYPIA INITIALIZE TTY PIA 
6F O1 CLR 1,X 
6F O3 CLR Sex 
86 01 LDA A #1 
A7 OO STA A X 
C6 O7 LDA B #7 
E7 O1 STAB 1,X 
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DISP 


00146 
00147 
00148 
00149 
00150 
00151 
00152 
00153 
00154 
00155 
00156 
00157 
00158 
00159 
00160 
00161 
00162 
00163 
00164 
00165 
00166 
00167 
00168 
00169 
00170 
00171 
00172 
00173 
00174 
00175 
00176 
00177 
00178 
00179 
00180 
00181 
00182 
00183 
00184 
00185 
00186 
00187 
00188 
00189 
00190 
00191 
00192 
00193 
00194 
00195 


3O5E 
3060 
3062 
3064 
3066 


3068 
306B 
3O6E 
3070 
3072 
3074 
3076 
3078 
307A 
307C 
SO7E 
3080 
3082 
3084 
3086 
3088 
308A 
308B 
308D 
308F 
3091 
3033 


3095 
3097 
309A 
309C 
3O9F 
3O0A2 
30A4 


30A7 
3OAA 


3OAC 
SOAF 
3OB 1 
30B4 
30B6 
30B9 
3OBB 
3OBE 
30CO 
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0O STA A X 
02 STA B 2,xX 

34 LDA A #$34 
03 STA A 3,X 
O02 STA A 2,X 
8020 LDX #$8020 INITIALIZE FFT CONTROL PIA’S 
8020 LDX #$8020 
04 LDA A #4 SET OUTPUTS TO $FF 
01 Sil Ame Amal ar BEFORE DEFINING DATA DIRECTION 
03 STARA 935 X 

FF LDA B #$FF 
Ooo STAB X 
O02 STA B 2,X 
01 CLR TeX 

iS CLR 3,X 

1C LDA B- #$1C DEFINE DATA DIRECTION 
OO STADB xX 
OE LDA B- #$0E 
02 Si Am Smee 

15 LDA B- #$15 SET CONTROL REGISTERS 
O1 SipAg BaealeX 
DEC B 

O03 STA A 3,€X 
A2 LDA A PHASE SET ADD/SUB/SWITCH CONTROL 
02 STA.A 2,X 

ie LDA A #$1C DETECT RISING EDGE FOR 1 HZ CLOCK 
23 STA A $23,X 
81 LDA A #$81 INIT. MC6840 CLIPPING COUNTER 
5009 STA A’ TIMER+1 
80 LDA B- #$80 

5008 STA B- TIMER 
5009 STA B- TIMER+1 
O03 LDA A #3 
5008 STA A TIMER 
BOOO LDX #DCO INIT SPECTRUM POINTER 
A4 STX DPTR 
CCCC LDX #$CCCC SCALE=204 .8 
9C STX SCALE 
CCO8 LDX #$CCO8 PRODUCES 10DB/VOLT ON CRT DISPLAY 
9E STX SCALE+2 
0000 LDX #$0000 

98 STX FACTOR FACTOR=3.0 
CcO02 LDX #$COO2 

9A STX FACTOR+2 
0000 LDX #$0000 OFFSET=0O 
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DISP 


00196 
00197 
00198 
00199 
00200 
00201 
00202 
00203 
00204 
00205 
00206 
00207 
00208 
00209 
00210 
00211 
00212 
00213 
00214 
00215 
00216 
00217 
00218 
00219 
00220 
00221 
00222 
00223 
00224 
00225 
00226 
00227 
00228 
00229 
00230 
00231 
00232 
00233 
00234 
00235 
00236 
00237 
00238 
00239 
00240 
00241 
00242 
00243 
00244 
00245 


30C3 


30C5 
30C8 
30CB 
30CD 
3OCF 
30D2 


3110 
3113 
3116 
3118 
311B 


311C 
311F 
3121 
3124 
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“AO STX OFFSET 

* 
311C LDX #IRQ SET ADDRESS OF INTERRUPT 
AQOO STX $A000 HANDLING ROUTINE 
cc LDA A #$CC SET HEADER FOR TAPE STORAGE 
5B STA A HEAD 
040D LDX #$040D 
5C STX HEAD+ 14 

bd 

* START FFT 

Bd 

START SEI 
OOB6 CLR DMAFLG 
ic LDA A #$1iC CLEAR FFT PROCESSOR 
8020 STA A’ FCNTRL 
10 LDA A #$10 ENABLE FFT PROCESSOR 
8020 STA A’ FCNTRL 
03 LDA A #3 INIT. PASS COUNTER 
57 STA A PASS 
A2 LDA A PHASE SET LOAD/ACCUMULATE TO LOAD 
08 ORA A #8 
8022 STA A PCNTRL 

* 
8020 PASS1 LDA A’ FCNTRL WAIT FOR FFT DONE FLAG 
8021 PASS2 LDA A FFTFLG f 
FB BPL PASS2 WAIT FOR 3 FFT’S 
0057 DEC PASS 
F3 BNE PASS1 
A2 LDA B- PHASE SET LOAD/ACCUM. TO ACCUM. 
3110 JSR CLK 
8020 LDA A’ FCNTRL 
0054 CLR COUNT 1 CLEAR INTEGRATION COUNTER 
0055 CLR COUNT 2 
O02 LDA A #2 
56 STA A COUNTS 

CLI ENABLE INTERRUPTS 

3523 JUMP DISP GO TO SPECTRUM DISPLAY ROUTINE 

* 
8022 CLK LDA A PCNTRL SUBROUTINE TO WAIT 
8023 CLK1 LDA A CLKFLG FOR FFT CLOCK TRANSITION 
FB BPL CLK1 
8022 STA B- PCNTRL 

RTS 

* 

* INTERRUPT HANDLING ROUTINE; 

* 
8021 IRQ TST FFTFLG CHECK FFT INTERRUPT 
24 BMI IRQFFT 
8010 LDA A ACIACR CHECK TAPE INTERRUPT 
oc BMI IRQA 
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00256 
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00258 
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00261 
00262 


3126 
3129 
312B 
3120 
312F 
3132 


3135 
3138 
313B 
3S13E 
3140 


3141 
3142 


3145 
3148 
314A 
314C 
314E 
314F 
3151 
3153 
3156 
3158 
315A 
315C 
315E 
3160 
3163 
3166 


3167 
3169 
316B 
316E 
316F 
3172 


3175 
3178 
317A 
317D 
3S17F 
3182 
3185 


3188 
318B 
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IRQA 
* 


IRQTMR 


* 


NEWIT 


ak 


IRQFFT 


NE 


RET 


DMA 


DMA 1 


LDA 
AND 
CMP 
BEQ 
JMP 
JUMP 


LDA 
LDA 
DEC 
BLE 
RTI 


LDA 


CLR 


LDX 
STX 
LDX 
STX 
LDX 
STX 
STX 


LDA 
ASL 


PPP 


b> p> 


HZ+1 
#$7F 
#$5C 
IRQTMR 
IRQHZ 
IRACIA 


TIMER+1 
TIMER+6 
TMRDEL 
NEWIT 


ITER1 


FCNTRL 
DMAFLG 
DMA 

COUNT 2 


COUNT 2 
NE 
COUNT 1 
INT+1 
RET 
COUNT 1 
INT 
RET 
DMAFLG 
CHFREQ 


PHASE 
#8 
CLK 


FCNTRL 
DMAFLG 


TIMER+2 
CLIPS+2 
TIMER+4 
CLIPS 
#$FFFF 
TIMER+2 
TIMER+4 


FFTFLG 


Page 6 


CHECK IF 1 HZ INTERRUPTS ENABLED 


ITERATION DELAY ROUTINE 
CLEAR INTERRUPT FLAG 


IF SUFFICIENT DELAY,BEGIN 


NEW DELETION ITERATION 


FFT INTERRUPT ROUTINE 


CLEAR INTERRUPT FLAG 
INCREMENT INTEGRATION COUNTER 


COMPARE TO REQUIRED COUNT 


IF EQUAL, START DMA 


SET LOAD/ACCUM. TO LOAD 


ENABLE DMA TRANSFER 


READ NUMBER OF CLIPS COUNTED 


WAIT FOR FFT DONE FLAG 
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31B1 


31B3 
31B5 
31B7 
31BA 
31BD 
31BE 
3 1BF 
31CO 
31C1 
31C4 
31C6 
31C9 
31CC 
31CE 
31CF 
31D0 


31D1 
31D4 
31D6 


31D8 
31DB 
31DD 
3 1DF 
31EO 
31E3 
31E5 
31E8 
31EB 


SiEE 
31F 1 
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FA BPL DMA 1 
8020 LDA A’ FCNTRL CLEAR INTERRUPT FLAG 
10 LDA A #$10 
8020 STA A FCNTRL DISABLE DMA TRANSFER 
A2 LDA B- PHASE SET LOAD/ACCUM. TO ACCUM. 
3110 JSR CLK 
0054 CLR COUNT 1 CLEAR INTEGRATION COUNTER 
0055 CLR COUNT2 
O02 LDA A #2 
56 STA A COUNT3 
CEI 
3C PUNCH LDA A #$3C 
8007 STA A’ RDRCTL START TAPE RECORDER MOTOR 
009 1 CLR CKSM INIT. CHECKSUM AND 
OO05B LDX #HEAD TAPE OUTPUT POINTER 
92 STX CASPTR 
* 
20 LDA A #$20 SET EXP=32 
A3 STA A EXP 
B800 LDX #CO CONVERT CO SPECTRUM FROM 
327F THREE JSR FPCONV INTEGER TO FLOATING POINT 
INX 
INX 
INX 
INX 
BAOO CPX #QUAD 
O03 BNE FP 1 
OOA3 DEC EXP SET EXP=31 TO DIVIDE BY 2 
COOO FP 1 CPX #QUAD+$600 CONVERT QUADRATURE AND AUTO 
EC BNE THREE SPECTRA TO FLOATING POINT 
NOP 
NOP 
NOP 
* 
33CA JSR ALLCON PRODUCE LOG DISPLAY VALUES 
TF LDA A #$7F SET OLDN=127 
A7 STA A OLDN 
o 
3FOO DINIT LDX #DVCTR1 INIT. DELETION VECTORS TO 1’S 
01 LDA A #1 3 
oO DINIT1 STA A xX 
INX 
3F80 CPX #DVCTR1+$80 
F8 BNE DINIT1 
OO6C CLR ITER CLEAR ITERATION COUNTER 
3604 JSR CLPOUT OUTPUT NO. OF CLIPS AND L.O. FREQ. 
3656 JSR FRQOUT TO CHART RECORDER 
* 
OO6C ITER1 INC ITER DELETION ITERATION ROUTINE 
3FOO LDX #DVCTR1 
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3234 
3237 
3239 
323C 
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3241 
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3250 
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3257 
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01 
0O 


FA 
SF7F 
01 
00 


CLRDV1 


CLRDV2 


CLR 
LDA 
STA 
INX 
DEC 
BNE 
LDX 
LDA 
STA 
DEX 
DEC 
BNE 


LDA 


‘STA 


LDX 
STX 
LDX 
STX 
LDX 
STX 
JSR 


STX 
JSR 
JSR 


LDX 
STX 
LDX 
STX 
LDX 
STX 
JSR 


STX 
JSR 
JSR 
JSR 


LDA 
CMP 
BLT 
BNE 
LDA 
STA 
LDA 
STA 


rPwP 


PrP 


PPP yp 


CLRDV1 
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Page 8 


#DVCTRI+$7F 


# 
Xx 


CLRDV2 


#3 
TIMER 


#CO 
DATPTR 
#MEAN1 
RESPTR 
#DVCTR1 
DELPTR 
MEAN 
#DVCTR2 
DELPTR 
COMP 
DVAND 


#QUAD 
DATPTR 
#MEAN2 
RESPTR 
#DVCTR1 
DELPTR 
MEAN 
#DVCTR2 
DELPTR 
COMP 
MNOUT 
DVAND 


ITER 
#3 
ITERS 
ITER2 
#$14 
H2+1 
#$3D 
ACIACR 


DISABLE DELAY TIMER 


INIT. POINTERS FOR CO SPECTRUM 


FIND MEAN AND STD. DEV. 

OF NON-DELETED CHANNELS 

DELETE CHANNELS WHICH DEVIATE 

FROM MEAN BY TOO MUCH 

COMBINE CO AND QUAD. DELETION VECTORS 


INIT. POINTERS FOR QUADRATURE 
SPECTRUM AND PERFORM DELETIONS 


OUTPUT MEANS TO CHART RECORDER 
COMBINE CO AND QUAD. DELETION VECTORS 


WAIT FOR 3 ITERATIONS BEFORE 
SENDING DATA TO TAPE RECORDER 


ENABLE DATA SENDING 
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ITER2 


ITERS 


ITER4 


ITERS 


* 
* 


FPCONV 


* 


OO4F NEG 


16 
0050 
14 
0051 
12 
00 
10 
01 


CMP 
BEQ 
STA 
LDA 
CMP 
BLT 
JSR 


RTI 


LDX 
STX 
LDA 
STA 
LDA 
STA 
RTI 


CLR 
STX 
LDA 
AND 
STA 
STX 
LDX 
LDA 
STA 
LDA 
LDX 
LDA 
STA 
STA 
LDA 
STA 
LDA 
STA 
LDA 
STA 
BPL 


NEG 
BCS 
NEG 


NEG 
BCS 
NEG 
BCS 
NEG 


PPPYp r>PPw 


Ppp 


PPP 


PPPrPrprroarwa 
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OLDN 
ITER4 
OLDN 
ITER 
ITMAX 
ITERS 
NOUT 
PRINT 


#SFFFF 
TIMER+6 
#$43 
TIMER 
#2 
TMRDEL 


SIGN 
LOW 
LOW 
#$F7 


HIGH 


WERE MORE CHANNELS DELETED? 

HAS MAX. NO. OF ITERATIONS 

BEEN REACHED? 

OUTPUT NO. OF CHANNELS TO RECORDER 
PRINT RESULTS ON TTY 


START ITERATION DELAY TIMER 


SUBROUTINE TO CONVERT 48-BIT 
INTEGER TO AM9511 32-BIT FLOATING 
POINT FORMAT 
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32F9 
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3301 
3303 
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330A 
330C 
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3310 
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3315 
3317 
3319 
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OE 
02 
Oc 
* 
0050 COM2 
0051 COM3 
00 COM4 
01 COMS 
02 COM6 
004A SGNST 
* 


A3 ZSTRT 


02 LMS 


03 SHIFT 


F3 SHIFT1 


03 DONE 


BCS 
NEG 
BRA 


COM 
COM 
COM 
COM 
COM 
INC 


LDA 
STA 
LDA 
LDA 
BNE 
LDA 
STA 
LDA 
STA 
LDA 
STA 
LDA 
STA 
LDA 


CLR 


LDA 
SUB 
STA 
DEC 
BNE 
CLR 
RTS 


DEC 
ASL 
ROL 
ROL 


BPL 
LDA 
AND 
TST 
BEQ 
ORA 
STA 
RTS 


Oprppp 


ODooaonannnaann 


roOwWw 


COM6 
SGNST 


LS2 


SHIFT 1 


LS3 


LS2 


Sit 
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* 
OO LDTOS LDA A X SUBROUTINE TO LOAD FLOATING 
5000 STA A TOS POINT NUMBER FROM LOCATION 
01 LDA A 1,X INDICATED BY X TO TOS 
5000 STA A TOS 
02 LDA A 2,xX 
5000 STA A TOS 
03 LDA A 3,xX 
5000 STA A TOS 
RTS 
* 
CL4TOS CLR A SUBROUTINE TO CLEAR 4 BYTES 
5000 STA A TOS ON TOS 
5000 STA A TOS 
5000 STA A TOS 
CL1TOS CLR A SUBROUTINE TO CLEAR 1 BYTE 
5000 STA A TOS ON TOS 
RTS 
* 
DISCON CLR B SUBROUTINE TO CONVERT FLOATING 
OO LDA A X POINT SPECTRUM VALUE TO 
O01 BEQ D1 LOGARITHMIC CRT DISPLAY VALUE 
INC B 
5000 D1 STA A TOS 
01 Le Ae Aen 
O1 BEQ D2 
INC B 
5000 D2 STA A TOS 
O02 ED Am Aga 2x 
O01 BEQ D3 
INC B 
5000 D3 STA A TOS 
03 LDA A 3,X 
O1 BEQ D4 
INC B 
5000 D4 STA A TOS 
TST B 
47 BEQ DEXIT 
OO4A CLR SIGN 
TST A 
08 BPL DPOS 
O04A INC SIGN 
15 LDA A #4$15 CHSF 
33A3 JSR CMD 
08 DPOS LDA A #8 LOG 
33A3 JSR CMD 
oosc LDX #SCALE 
331C JSR LDTOS 
12 LDA A #$12 FMUL 
33A3 JSR CMD 
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DISP 


00546 
00547 
00548 
00549 
00550 
0055 1 
00552 
00553 
00554 
00555 
00556 
00557 
00558 
00559 
00560 
0056 1 
00562 
00563 
00564 
00565 
00566 
00567 
00568 
00569 
00570 
0057 1 
00572 
00573 
00574 
00575 
00576 
00577 
00578 
00579 
00580 
0058 1 
00582 
00583 
00584 
00585 
00586 
00587 
00588 
00589 
005390 
00591 
00592 
00593 
00594 
00535 


3382 
3384 
3387 
3389 
338C 
338E 
3391 
3393 
3396 
3397 
3399 
339C 
339F 
33A1 
33A3 
33A6 
33A9 
33AB 


33AC 
33AF 
33B1 
33B4 
33B5 
33B7 
33BA 
33BB 
33BD 
33BE 
33BF 


33C 1 
33C3 
33C4 
33C5 
33C6 
33C7 
33C9 


33CA 
33CD 
33CF 
33D2 
33D4 
33D7 
33DA 
33DC 
33DF 
33E 1 
33E4 
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3340 


5000 


0O 
5000 


01 


48 


52 


D5 


CMD 


BUSY 


DEXIT 


CHCON 


* 
INCHAN 


* 


ALLCON 


Al 


A2 


LDA 


LDX 


RTS 


LDX 
STX 
LDX 
STX 
JSR 
CPX 
BNE 
JSR 
LDA 
JSR 
JSR 


PPPpPpp > 


> 


PPP 


PPPrYprpp 


#$1F 
CMD 
OFFSET+1 
TOS 
OFFSET 
TOS 
#$6D 
CMD 


DISCON 
DISPTR 
TOS 


X 
TOS 


1,X 


DISPTR 


CHAN 


CHAN 


#$BO0O 
DISPTR 
#$B800 
CHAN 
CHCON 
#$BCOO 
Al 
LOTOS 
#$15 
CMD 
STTOS 
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FIXS 


SSUB 


CHSS 


SUBROUTINE TO CONVERT A SPECTRAL 
COMPONENT TO LOG DISPLAY VALUE 


SUBROUTINE TO CONVERT ALL 4 
SPECTRA} TO LOG DISPLAY VALUES 


CHSF 
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DISP 


00596 
00597 
00598 
005399 
00600 
00601 
00602 
00603 
00604 
00605 
00606 
00607 
00608 
00603 
006 10 
O06 11 
006 12 
00613 
O06 14 
006 15 
006 16 
006 17 
006 18 
006 19 
00620 
00621 
00622 
00623 
00624 
00625 
00626 
00627 
00628 
00629 
00630 
00631 
00632 
00633 
00634 
00635 
00636 
00637 
00638 
00639 
00640 
0064 1 
00642 
00643 
00644 
00645 


33E7 
S3EA 
33ED 
S3EF 


33FO 
33F2 
33F4 
33F6 
33F8 
33FB 
S3FE 


3401 
3403 
3405 
3407 
3409 
340C 
340E 
3411 
3413 
3416 
3419 
341B 
341E 
3420 
3423 
3425 
3428 
342A 
342D 


3430 
3433 
3435 
3436 
3438 
343A 
343C 
343E 
3441 
3444 
3446 
3449 
344B 


344E 
3450 
3453 
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33AC JSR CHCON 
COOO CPX #$COOO 
ED BNE A2 
RTS 

* 5 
44 MEAN LDX DATPTR SUBROUTINE TO FIND MEAN AND 
52 STX CHAN STD. DEV. OF SPECTRAL COMPONENTS 
46 LDX DELPTR FOR WHICH DELETION VECTOR=1 
58 STX DEL 
OO5A CLR NUM DATPTR SPECIFIES BEGINNING OF SPECTRUM 
3331 JSR CL4TOS 
3331 JSR CL4TOS DELPTR SPECIFIES BEGINNING OF 

* DELETION VECTOR 
58 M1 LDX DEL 
00 LDA A xX RESPTR SPECIFIES WHERE TO STORE RESULTS 
29 BEQ M2 
52 LDX CHAN 
331C JSR LDTOS 
10 LDA A #$10 FADD 
33A3 JSR CMD 
19 LDA A #$19 X CHF 
33A3 JSR CMD 
331C JSR LDTOS 
17 LDA A #$17 PTOF 
33A3 JSR CMD 
12 LDA A #$12 FMUL 
33A3 JSR CMD 
10 LDA A #$10 FADD 
33A3 JSR CMD 
19 LDA A #$19 XCHF 
33A3 JSR CMD 
OOSA INC NUM 

* 
33C1 M2 JSR INCHAN 
59 LDA A DEL+1 

INC A 

59 STA A DEL+1 
7F AND A #$7F 
cs BNE M1 
SA LDA A NUM 
5000 STA A TOS 
333B JSR CL1TOS 
1D LDA A #$1D ARIES 
33A3 JSR CMD 
13 LDA A #$13 FDIV 
33A3 JSR CMD 

* 
17 LDA A #$17 PTOF 
33A3 JSR CMD 
42 LDX RESPTR 
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3462 
3464 
3467 
346A 
346C 
S46F 
3471 
3474 
3476 
3479 
347B 
347D 
S47E 
3481 
3484 
3486 
3489 
348B 
348E 
3430 
3493 
3494 
3495 
3496 
3497 
3498 
349B 
349D 
34A0 
34A2 
34A5 
34A7 
S4AA 
34AC 


34AD 
S4AF 
34B 1 
34B3 
34B5 
34B7 
34B8 
34B9 
34BA 
34BB 


Motorola M68SAM Cross-Assembler 


STTOS 


COMP 


rPPPpPp > 


> 


PPrPPrPrPprKPrYpp 


STTOS 
#$17 
CMD 
#$12 
CMD 
NUM 
TOS 
CL1TOS 
#$1D 
CMD 
#$12 
CMD 
#$11 
CMD 
NUM 
4,X 


TOS 
CL1iTOS 
#$1D 
CMD 
#$13 
CMD 
#$01 
CMD 


TOS 
3,X 
TOS 
2,X 
TOS 
1,X 
TOS 


DATPTR 
CHAN 
DELPTR 
DEL 
RESPTR 
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PTOF 


FMUL 


FLTS 
FMUL 


FSUB 


Bes 
FDIV 


SQRT 


SUBROUTINE TO STORE TOS AT 
LOCATIONS INDICATED BY X 


SUBROUTINE TO COMPARE CHANNELS 
TO MEAN AND DELETE THOSE BEYOND 
FACTOR*SDEV 
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DISP 


00696 
00697 
00698 
00699 
00700 
00701 
00702 
00703 
00704 
00705 
00706 
00707 
00708 
00703 
00710 
007 11 
00712 
00713 
007 14 
00715 
007 16 
007 17 
007 18 
00719 
00720 
00721 
00722 
00723 
00724 
00725 
00726 
00727 
00728 
00729 
00730 
00731 
00732 
00733 
00734 
00735 
00736 
00737 
00738 
00739 
00740 
00741 
00742 
00743 
00744 
00745 


34BC 
34BF 
34C2 
34C5 
34C7 
34CA 
34CD 


34D0 
34D3 
34D6 
34D8 
34DB 
34D0D 
34E0 
3S4E2 
34E5 
S4E6 
34E8 
S4EA 
34ED 
34EF 
34F2 
34F4 
B4F6 
34F7 
34F9 
34FB 
34FC 
S4FE 
3501 
3503 
3505 
3507 


3509 
350B 
350D 
350E 
3510 
3511 
3513 
3515 
3517 
3518 
351A 
351C 
351E 
3520 
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331C JSR LDTOS 
0098 LDX #FACTOR 
331C JSR LDTOS 
12 LDA A #$12 FMUL 
33A3 JSR CMD 
0094 LDX #MAX 
3498 JSR STTOS 
* 
0094 Ci LDX #AMAX 
331C JSR LDTOS fs 
42 LDX RESPTR 
331C JSR LDTOS 
52 LDX CHAN 
331C JSR LDTOS 
11 LDA A #$11 FSUB 
33A3 JSR CMD 
ASL A 
O05 BPL C2 
15 LDA A #$15 CHSF 
33A3 JSR CMD 
iit C2 LDA A #$11 FSUB 
33A3 JSR CMD 
58 LDX DEL 
0O CLR Xx 
ASL A 
02 BMI C3 
OO INC X 
C3 INX 
58 STX DEL 
33C1 JSR INCHAN 
38) LDA A DEL+1 
tr AND A #$7F 
00 CMP A #$00 
C7 BNE C1 
* 
* 
46 COMP2 LDX DELPTR DELETION OF ADJACENT CHANNELS 
oO1 LDA B 1,X 
c4 TBA 
O02 LDA B 2,X 
TST A 
04 BNE cs 
0O CLR X 
02 CLR 2,X 
c5 INX 
58 STX DEL 
59 LDA A DEL+1 
Ue AND A #$7F 
7E CMP A #$7TE 
EB BNE C4 
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DISP 


00746 
00747 
00748 
00749 
00750 
00751 
00752 
00753 
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00755 
00756 
00757 
00758 
00759 
00760 
0076 1 
00762 
00763 
00764 
00765 
00766 
00767 
00768 
00769 
00770 
0077 1 
00772 
00773 
00774 
00775 
00776 
00777 
00778 
00779 
00780 
0078 1 
00782 
00783 
00784 
00785 
00786 
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00730 
0079 1 
00792 
00793 
00794 
00795 


3522 


3523 
3525 
3527 
3529 
352B 
352D 
352F 
3532 
3534 
3536 
3538 
353B 
353D 
353F 
3541 
3543 
3546 
3548 
354A 
354C 
354E 
3551 
3553 
3555 
3557 
3559 
355C 
SS5E 
3560 
3562 
3564 
3566 
3569 
356C 
356F 
3572 


3575 
3577 
3579 
357B 
357C 
S57E 
3581 
3583 
3586 


Motorola M68SAM Cross-Assembler Page 16 
RTS 

* 

* 
50 DISP BSR SWEEP CRT DISPLAY ROUTINE 
FC BMI DISP 
4C DBNC BSR SWEEP READS KEYBOARD, DEBOUNCES KEYS 
FC BPL DBNC AND INTERPRETS KEYS 
48 BSR SWEEP 
F8 BPL DBNC 
8008 KEY LDA A KEYBD 
1F AND A #$1F 
18 STRNG1 CMP A ¥#$18 ONE - DISPLAY CO SPECTRUM 
O7 BNE STRNG2 
BOOO LDX #DCO 
A4 STX DPTR 
E4 BRA DISP 
19 STRNG2 CMP A #$19 TWO - DISPLAY QUADRATURE SPECTRUM 
O7 BNE STRNG3 
B100 LDX #DQUAD 
A4 STX DPTR 
DS BRA DISP 
1A STRNG3 CMP A #$1A THREE - DISPLAY AUTO1 
O7 BNE STRNG4 
B200 LDX #DAUTO1 
A4 STX DPTR 
GE BRA DISP 
1B STRNG4 CMP A #$1B FOUR - DISPLAY AUTO2 
O7 BNE PTEST 
B300 LDXx #DAUTO2 
A4 STX DPTR 
C3 BRA DISP 
17 PTEST CMP A #$17 BLANK - STOP AND GO TO MONITOR 
BF BNE DISP 
18 LDA A #$18 
8020 STA A FCNTRL STOP FFT PROCESSOR 
30D4 LDX #START 
A048 STX $A048 SET RESTART POINTER 
AO7TF LDS #$AOTF SET STACK POINTER 
E11C JMP $E11C GO TO MONITOR 

* 

* 
A4 SWEEP LDX DPTR SUBROUTINE TO GENERATE 
O7 LDA A #7 X AND Y AXES FOR CRT 
AG STA A RAMPH 

RMPH CLR B 
oO YDATA LDA A X 
D204 STA A DACS3 DAC3 - Y AXIS SPECTRUM 
O01 LDA A 1,X 
D205 STA A DACS3+1 

INX 
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DISP 


00796 
00797 
00798 
00739 
00800 
00801 
00802 
00803 
00804 
00805 
00806 
00807 
00808 
00809 
008 10 
008 11 
008 12 
008 13 
008 14 
00815 
008 16 
008 17 
008 18 
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00822 
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00825 
00826 
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00844 
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3587 
3588 
358A 
358D 
358F 
3592 
3593 
3595 
3597 
3598 
359A 
359D 
359F 
35A 1 
35A3 
35A5 
35A8 
SSAA 
35AD 
35BO 


35B3 
35B6 
35B9 


35BA 
35BB 
35BE 
35CO 
35C2 
35C4 
35C5 
35C7 
35C8 
35CB 
35CD 


35CE 
35D1 
35D4 
35D7 
35DA 
35DD 
35E0 
35E3 
35E6 
35E7 
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35EC 
35ED 
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INX 
AG LDA A RAMPH 
D202 STA A DAC2 DAC2 - X AXIS RAMP 
04 RMPL SUB B #4 
D203 STA B DAC2+14 
TBA 
1F AND A #$1F 
F6 BNE RMPL 
TST B 
E2 BNE YDATA 
OOA6 DEC RAMPH 
AS LDA A RAMPH 
F7 CMP A #$F7 
D8 BNE RMPH 
O7 LDA A #7 
D202 STA A DAC2 
EE LDA A AFF 
D203 STA A DAC2+1 
D204 STA A DAC3 
D205 STA A DAC3+1 
* 
8008 PRESS LDA A KEYBD CHECK KEYBOARD 
8009 LDA A KEYBD+1 
RTS 
ok 
DVAND CLR B SUBROUTINE TO COMBINE 
3FOO LDX ADVCTR1I CO AND QUAD DELETION VECTORS 
00 DVAND1 LDA A xX 
80 AND A $80,X — 
O01 BEQ DVAND2 
INC B ; 
00 DVAND2 STA A X 
INX 
3F80 CPX #DVCTR1+$80 
Fi BNE DVAND 1 
RTS 
* 
OO6D MNOUT LDX #MEAN1 SUBROUTINE TO OUTPUT CO AND QUAD 
3340 JSR DISCON MEANS TO CHART RECORDER 
D206 LDX #DAC4 
35E3 JSR DACOUT 
0076 LDX AMEAN2 
3340 JSR DISCON 
D208 LDX #DACS 
5000 DACOUT LDA A_ TOS DAC OUTPUT SUBROUTINE 
COM A 
00 STA A X 
5000 EDARA LOS 
COM A 
01 STA A 1,X 
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DISP 


00885 


00886 
00887 
00888 
00889 


SSEF 


35FO 
35F2 
35F4 
35F5 
35F6 
35F7 
35F8 
35F9 
35FA 
35FB 
35FC 
35FD 
3600 
3603 


3604 
3607 
3609 
360B 
360E 
3610 
3613 
3615 
3618 
361B 
S61E 
3620 
3623 
3626 
3629 
362B 
362E 
3630 
3633 
3636 


3639 
363A 
363B 
363C 
363D 
363E 
363F 
3640 


3641 
3643 
3646 
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A7 
40 


03 
5000 
O02 


* 


NOUT 


CLPOUT 


KMUL 


KSUB 


* 
LR4TOS 


RTS 


LDA 


RTS 


LDA 
STA 
LDA 


OrPraroaoarwrwrarp 


PPP 


OLDN 
#$40 


DAC6 
DAC6+1 


#CLIPS 
LR4TOS 
#$34 
CMD 
#$1C 
CMD 

#8 

CMD 
#KMUL 
LDTOS 
#$12 
CMD 
#KSUB 
LDTOS 
#$11 
CMD 
#$1F 
CMD 
#DAC7 
DACOUT 
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SUBROUTINE TO OUTPUT NO. OF 
CHANNELS NOT DELETED TO RECORDER 


SUBROUTINE TO OUTPUT LOG OF 
NO. OF CLIPS TO RECORDER 
CHSD 

FLTD 


LOG 


FMUL 


FSUB 


FIXS 


0,0,$80,$0A 


0,0,$80,$0C 


SUBROUTINE TO LOAD TOS FROM 
LOCATIONS INDICATED BY X 
IN REVERSE ORDER 
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009 11 
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009313 
009 14 
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009 16 
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00919 
00920 
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00822 
00923 
00924 
00325 
00926 
00927 
00328 
00929 
00930 
00931 
00332 
00933 
00334 
00935 
00936 
00937 
00938 
00339 


3648 
364B 
364D 
3650 
3652 
3655 


3656 
3659 
365B 
365D 
S65F 
3661 
3663 
3665 
3667 
366A 
366C 
366E 
3671 
3674 
3676 
3678 
367B 
S67E 


3681 
3683 


3685 
3686 
3688 
3689 
368B 


368D 
368E 
368F 
36390 
3691 
3693 
3696 
3699 
369B 
369E 
36A1 
36A4 
36A6 
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05 


OF 
04 


OA 
5000 
333B 


33A3 
5000 
333B 


33A3 


LR2TOS 
LR1ITOS 


* 


FRQOUT 


x 


FRQSUB 
FRQMUL 
* 


UNPCK 


* 
* 


DIGITH 


DIGITL 


STA 
LDA 
STA 
LDA 


RTS 


PSH 
BSR 
PUL 
AND 
BRA 


LSR 


LSR 
LSR 
LDA 
STA 
JSR 
LDA 
JSR 
STA 
JSR 
LDA 
JUMP 
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TOS 
1,X 
TOS 
X 

TOS 


CL4TOS 
FREQ 
#$3F 
UNPCK 
FREQ+1 
UNPCK 
FREQ+2 
DIGITH 
#FRQSUB 
LR2TOS 
#$6D 
CMD 
#FROMUL 
LR2TOS 
#$6E 
CMD 
#DAC8 
DACOUT 


$0000 
$0001 


DIGITH 


#$OF 
DIGITL 


#$OA 
TOS 
CLiTOS 
#$6E 
CMD 
TOS 
CLiITOS 
#$6C 
CMD 


SUBROUTINE TO OUTPUT L.O. 
FREQUENCY TO RECORDER 


SSUB 


SMUL 


SUBROUTINE TO CONVERT PACKED 
BCD NUMBER TO BINARY 


SMUL 


SADD 
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00959 
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0096 1 
00962 
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00966 
00967 
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00969 
00970 
0097 1 
00972 
00973 
00974 
00975 
00976 
00977 
00978 
00979 
00980 
0098 1 
00982 
00983 
00984 
00985 
00986 
00987 
00988 
oo0s8s9g 


36A9 
36AB 
36AD 
S6AF 
36B2 
36B4 
36B6 
36B9 
36BC 


3S6BE 
36C1 
36C3 
36C6 
36C8 
36CA 
36CB 
36CD 
36CE 


S6CF 
36D1 
36D4 
36D6 
36D8 
36D9 
36DB 


36DC 
36DE 
36E0O 
36E3 
36E4 
36E6 


36E7 
36E9 
36EC 
SG6EE 
S6F 1 


36F2 
36F5 
36F7 
36F9 
36FA 
36FC 
S6FE 
3700 
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92 


0O 
8011 


91 


92 


91 


8011 


IRACIA 


SPECT 


WRITE 


CSOUT 


* 


CASSTP 


TRQHZ 


LDX 
LDA 
BNE 
CPX 
BNE 
LDA 
STA 
LDX 
STX 


CPX 
BEQ 
CPX 
BEQ 
LDA 
LSR 
BCS 
INX 
INX 


LDA 
STA 
ADD 
STA 
INX 
STX 
RTI 


LDA 
SUB 
STA 
INX 
STX 
RTI 


LDA 
STA 
LDA 


RTI 


LDA 
LDA 
ADD 
DAA 
STA 
CMP 
BCS 
CLR 
LDA 


nwo nnowowa 


nonnnw 


bP yp 


Cross-Assembler Page 20 


CASPTR 
CASPTR 
SPECT 
#HEAD+17 
WRITE 
#$1C 
HZ+1 

#CO 
CASPTR 


#$COOO 
CSOUT 
#$CO0O1 
CASSTP 
CASPTR+1 


WRITE 


X 
ACIATR 
CKSM 
CKSM 


CASPTR 


CKSM 
#$DD 
ACIATR 


CASPTR 


#$1D 
ACIACR 
#$34 
RDRCTL 


HZ 
TIME+2 
#1 


TIME+2 
#$60 

IRQHZ 1 
TIME+2 
TIME+1 


ROUTINE FOR TAPE RECORDER OUTPUT 
CHECK IF HEADER OR SPECTRA 
ENABLE 141 HZ INTERRUPTS 

WHEN HEADER COMPLETED 

SET POINTER TO CO SPECTRUM 


CHECK IF ALL SPECTRA COMPLETED 


CHECK IF ADDRESS IS EVEN OR ODD 


IF EVEN, INC. X TWICE TO SKIP 
LAST TWO BYTES OF SPECTRAL COMPONENT 


OUTPUT BYTE INDICATED BY X 


OUTPUT CHECKSUM 


DISABLE TAPE INTERRUPTS 
AND TURN TAPE OFF 


1 HZ CLOCK ROUTINE 
COUNTS SECONDS, MINUTES, AND HOURS 
FOR TIME OF DAY 
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DISP 
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3705 
3707 
3708 
370A 
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3711 
3713 
3715 
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3718 
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371C 
37 1F 


3720 
3721 
3722 
3723 
3724 
3725 
3726 
3727 
3728 
3729 
372A 
372B 
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372D 
S72E 
372F 


3730 
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3740 
3743 
3745 
3748 
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374D 


374F 
3752 
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3757 
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O01 ADD #1 
DAA 
60 STA TIME+1 
60 CMP #$60 
11 BCS IRQHZ 1 
0060 CLR TIME+1 
SF LDA TIME 
O01 ADD #1 
DAA 

SF STA TIME 
24 CMP #$24 
O03 BCS IRQHZ 1 
OO5F CLR TIME 

IRQHZ1 RTI 

* 

K10E8 FCB 5,$F5,$E1,0 

K10EG FCB 0O,$F,$42,$40 

K10E4 FCB 0,0,$27,$10 

K10E2 FCB 0,0,0,$64 

* 
3641 BYTE4 JSR LR4TOS SUBROUTINE TO PRINT NO. OF CLIPS 
34 LDA #$34 CHSD 
33A3 JSR CMD CONVERTS COMPLEMENT OF 4 BYTE 
3720 LDX #K10E8 INTEGER TO DECIMAL FOR PRINTING 
eh 7/ BYTE41 LDA #$37 PTOD 
33A3 JSR CMD 
3641 JSR LR4TOS 
Ns LDA #$2QF DDIV 
33A3 JSR CMD 
37 LDA #$37 PTOD 
33A3 JSR CMD 
16 BSR STKLST 

ok 
3641 JSR LR4TOS 
2E LDA #$2E DMUL 
33A3 JSR CMD 
2D LDA #$2D DSUB 
33A3 JSR CMD 
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375F 
3760 
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3778 
377A 
377D 
S77F 
3782 
3784 
3787 
3789 
378C 
378E 
3791 
3793 
3796 
3798 


379B 
379D 
37A0 
37A2 
37A5 
37A7 
S7AA 
37AC 
37AF 
37B2 
37B5 
37B7 
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37BC 
37BD 
37BE 
37CO 
37C2 
37C4 
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INX 
INX 
INX 
INX 
3730 CPX #K10E2+4 
D6 BNE BYTE41 
5000 STKLST LDA B_ TOS 
5000 LDA B- TOS 
5000 LDA B- TOS 
5000 LDA B TOS 
56 BRA INT2 
* 
OD PRINT LDA A #$0D SUBROUTINE TO PRINT TIME, 
37FO JSR OUTCH L.O. FREQUENCY, NO. OF CLIPS, 
OA LDA A #$OA AND NO. OF CHANNELS NOT DELETED 
37FO JSR OUTCH 
5F LDA A TIME 
37DC JSR OUT2H 
3A LDA A #$3A 
37FO JSR OUTCH 
60 LDA A TIME+1 
37DC JSR OUT 2H 
3A LDA A #$3A 
37FO JSR OUTCH 
61 LDA A TIME+2 
37DC JSR OUT2H 
20 LDA A #$20 SP 
37FO JSR OUTCH 
* 
65 LDA A FREQ 
37DC JSR OUT 2H 
66 LDA A FREQ+1 
37DC JSR OUT 2H 
67 LDA A FREQ+2 
37E2 JSR OUTHL 
20 LDA A #$20 SP 
37FO JSR OUTCH 
0068 LDX #CLIPS 
3730 JSR BYTE4 
20 LDA A #%$20 SP 
37FO JSR OUTCH 
A7 LDA B- OLDN 
* 
INTOUT CLR A SUBROUTINE FOR 1 BYTE BINARY 
INT 1 INC A TO DECIMAL CONVERSION AND PRINTING 
64 SUB B #$64 
FB BCC INT 1 
64 ADD B- #$64 
25 ADD A #$2F 
37FO JSR OUTCH 
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DISP 
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3812 
3814 
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F7 
02 


02 


AOOF 


INT2 
INTS 


OUT 2H 


OUTHL 


OUTHR 


OUTCH 


IOUT 


IOUT 1 


IOUT2 


Ios 


DEL1 


DE 


PSH 


PSH 


RTS 


isi) 
BPL 
INC 
DEC 
RTS 


OYrp 


> 


PYPPprprpryppTp > 
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#$OA 
INTS 
#$OA 
#A$QF 
OUTCH 


#$30 
OUTCH 


OUTHL 


OUTHR 


#$OF 
#$30 
#$39 
OUTCH 
#7 


XTMP 
ATTYPIA 
#$OA 

X 

DE 

DEL1 

X 


IOUT 1 
2,X 


10s 
DEL1 
XTMP 


OUTPUT HEX BYTE AS 2 ASCII 
CHARACTERS 


OUTPUT ONE ASCII CHARACTER 
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DISP 


01128 
01129 
01130 
01131 
01132 
01133 
01134 
01135 
01136 
01137 
01138 
01139 
01140 
01141 
01142 
01143 
01144 
01145 
01146 
01147 
01148 
01149 
01150 
01151 
01152 
01153 
01154 
01155 
01156 
01157 
01158 
01159 
01160 
01161 
01162 
01163 
01164 
01165 
01166 


381B 
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CHFREQ 


oe 
INFREQ 


SYNTH 


LDA 
STA 
ADD 
DAA 
STA 
LDA 
STA 
ADC 
DAA 
STA 
LDA 
STA 
ADC 
DAA 
STA 


CMP 
BNE 
LDX 
CPX 
BNE 


LDX 
STX 
LDA 
STA 


LDX 
STX 
LDA 
STA 
ORA 
STA 
AND 
STA 
RTS 


END 


rPPrpp PPP DYD rPPp 
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PPrPprYrrp 
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NFREQ+2 
FREQ+2 
FINC+2 


NFREQ+2 
NFREQ+1 
FREQ+1 
FINC+1 


NFREQ+1 
NFREQ 
FREQ 
FINC 


NFREQ 


FMAX 
SYNTH 
NFREQ+1 
FMAX+ 1 
SYNTH 


FMIN+1 
NFREQ+1 
FMIN 
NFREQ 


NFREQ+1 
LATCH+ 1 
NFREQ 
LATCH 
#$80 
LATCH 
AGSTF 
LATCH 
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Appendix 4 - Robust Estimation Program 


The following is a listing of the FORTRAN program used for 
robust maximum-1likelihood estimation on the recorded 


spectra. 
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PROGRAM READS RAW SPECTRA FROM UNIT 7, 

WRITES RESULTS ON UNIT 6. 

PERFORMS 3 SIGMA ROBUST ESTIMATION ON AUTO SPECTRA. 
ESTIMATES VARIANCE OF CROSS SPECTRA FROM FIRST DECILES 

OF AUTO SPECTRA. 

FINDS ROBUST LOCATION ESTIMATE FOR CROSS SPECTRA 

USING PROCEDURE COMBINING REJECTION, HUBER ESTIMATION 

AND BIWEIGHT ESTIMATION. 

RECORDS ESTIMATES PLUS NUMBERS OF OUTLIERS FOUND AND TOTAL 
ABSOLUTE POWER DELETED. 


VARIABLES USED: 

AMEAN(4) - MEANS 

CLIPS - COMPLEMENT OF NUMBER OF CLIPS RECORDED DURING SPECTRUM 

CLP - NUMBER OF CLIPS 

ERR - INDICATOR FOR RECORDING ERROR 

FILE - FILE NUMBER FOR SPECTRUM 

FREQ - SECOND LOCAL OSCILLATOR FREQUENCY 

ICRIT(4) - CRITERIA SATISFIED DURING ESTIMATIONS ITERATIONS 

IER - ERROR INDICATOR FOR SUBROUTINE BEIUGR 

IOPT(S) - CONTROL CODES FOR SUBROUTINE BEIUGR 

IPERM(108,2) - VECTORS FOR PERMUTATIONS OF CROSS SPECTRA 

DURING ORDERING 

NDEL(2) - NUMBER OF DELETION ITERATIONS 

NOUT(2) - NUMBER OF OUTLIERS 

PSI(600,2) - PSI FUNCTIONS FOR ROBUST MAXIMUM LIKELIHOOD ESTIMATION 
1 - BIWEIGHT 2 - HUBER 

R(512) - ARRAY FOR FORMAT CONVERSION OF SPECTRA 

S(108) - TEMPORARY STORAGE FOR ORDERED AUTO SPECTRA 

SDEV(4) - STANDARD DEVIATIONS 

SPECT(126,4) - RAW SPECTRA 1-IN PHASE, 2-QUADRATURE, 3,4-AUTO 

SPORD(108,2 ) - ORDERED CROSS SPECTRA 1-IN PHASE, 2-QUADRATURE 

SQ(108,2) - SQUARES OF AUTO SPECTRAL COMPONENTS 

STAT(5) - MEAN,MAX,MIN,MEDIAN, VARIANCE 

STRES(6,4) - ARRAY OF STATISTICS OF SPECTRA 

SUM(2) - SUMS FOR STANDARD DEVIATION CALCULATIONS 

THETA(2) - ROBUST LOCATION ESTIMATES 

TOLD(2) - PREVIOUS LOCATION ESTIMATES 

X(600) - INCREMENTAL X VALUES FOR PSI FUNCTION GENERATION 

TIME - TIME WHEN SPECTRUM WAS RECORDED 

DER(600,2) - DERIVATIVES OF PSI FUNCTIONS 

IDV(108) - VECTOR INDICATING DELETED POINTS 

REAL AMEAN(4),SDEV(4),X(600) , TOLD(2) ,SUM(2) 

INTEGER IPERM(108,2),ICRIT(4),NOUT(2),NDEL(2) 

INTEGER CLP,IOPT(S5)/1,0,0,1,1/,IER 

REAL $(108),STAT(5),STRES(6,4),SQ(108,2) 

INTEGER ERR,FILE,TIME,INT,FREQ,CLIPS 

REAL SPECT(126,4),R(512) 

EQUIVALENCE (ERR,R(1)), (FILE,R(2)),(TIME,R(3)),(INT,R(4)) 

EQUIVALENCE (FREQ,R(5)),(CLIPS,R(6)), (SPECT(1,1),R(7)) 

REAL PSI(600,2),DER(600, 2), THETA(2),SPORD( 108, 2) 

INTEGER IDV(108) 

COMMON PSI,DER, THETA, SPORD, IDV 


GENERATE HUBER AND BIWEIGHT INFLUENCE CURVES 
AND THEIR DERIVATIVES 


H=O. 1000E-01 
DO 5 I=1,600 
AOb = Cte) FH 
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DO 7 I=1,500 
PSI(I,1)=X(1I)*(1.0-(X(1)/5. )**2)**2 
DO 8 I=501,600 

PSI(I,1)=0.0 

DO 9 I=1,150 

PSI(I,2)=X(I) 

DO 10 I=151,600 

PSI(I,2)=1.5 

‘CALL DERIV(PSI(1,1),DER(1,1)) 

CALL DERIV(PSI(1,2),DER(1,2)) 


READ A RECORD 


READ(7,30, END=99)(R(I),1=1,512) 

FORMAT (32A4) 

IF(ERR.NE.O)GO TO 20 

CEP =—CRiPS =a 
WRITE(6,40)FILE,TIME,INT,FREQ,CLP 
FORMAT (“O07 (225 2X,26,2X,2Z6,2X,26.2X, 110) 


FIND MEAN,STD DEV,MEDIAN,FIRST DECILE, 
MAX AND MIN FOR AUTO SPECTRA 


DO 80 I1=3,4 
DO 50 J=1, 108 

S(JU)=SPECT(uJ+9,I1) 

CALL VSRTA(S, 108) 

IMSL LIBRARY SUBROUTINE FOR ORDERING ACCORDING TO SIZE 
CALL BEIUGR(S,108,IOPT,STAT,IER) 

IMSL LIBRARY SUBROUTINE FOR FINDING STATISTICS 
STRES(1,1)=STAT(1) 

STRES(2,1)=SQRT(STAT(5)) 

STRES(3,1)=STAT(4) 

STRES(4,1)=S(11) 

STRES(5,1)=S( 108) 

STRES(6,1)=S(1) 

WRITE(6,70)(STRES(UJ,I1I),U=1,6) 

FORMAT(’ ‘’,6(E10.4,2X)) 

CONTINUE 


CHANNEL DELETIONS FOR AUTO SPECTRA 


DO 170 I=1,2 

K=1I+2 

DO 90 J=1, 108 
SQ(JU,1I)=SPECT(U+9,K)**2 


INITIALIZE MEAN AND STD DEV 


AMEAN(K)=STRES(3,K) 
SDEV(K)=STRES(2,K) 
ITER=O 

NLAST=108 


DELETE CHANNELS 


DMAX=3. *SDEV(K) 
ITER=ITER+1 
DON1105U=1, 108 
IDV(U)=1 
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IF (ABS(SPECT(JU+9,K)-AMEAN(K) ).GT.DMAX)IDV(J) =O 
DELETE ADJACENT CHANNELS 


NOW=IDV(2) 
IF(IDV(1).EQ.0)IDV(2)=0 
DOM190 8 =25107, 
NEXT=IDV(J+1) 
IF(NOW.NE.O)GO TO 120 
IDV(JU-1)=0 

IDV(JU+1)=0 

NOW=NEXT 

CONTINUE 

IF (NOW.EQ.0)IDV(107)=0 


FIND MEAN,STD DEV,AND NO. OF CHANNELS AFTER DELETIONS 


NCHAN=0O 
SUMM=O. 

SUMS=O. 

DO 140 J=1,108 
IF(IDV(JU).EQ.0)GO TO 140 
SUMM=SUMM+SPECT (uJU+9,K) 
SUMS=SUMS+SQ(u,I) 
NCHAN=NCHAN+ 4 

CONTINUE 
AMEAN(K ) =SUMM/NCHAN 
SDEV(K)=SQRT( (SUMS- (SUMM* *2 ) /NCHAN ) /NCHAN ) 
IF(NCHAN.EQ.NLAST)GO TO 150 
IF(ITER.GT.20)GO TO 150 
NLAST=NCHAN 

GO TO 100 


WRITE RESULTS 


WRITE(6, 160)AMEAN(K),SDEV(K) ,NCHAN, ITER 
FORMAT(’ ’,£10.4,2X,E£10.4,2X,13,2x,12) 
CONT INUE 


FIND MEAN, STD DEV, MEDIAN,MAX AND MIN FOR CROSS SPECTRA 


DOS21OFT=172 
DO 200 J=1, 108 

SPORD(JU,1)=SPECT(uU+9,1) 

IPERM(J,I)=uJ 

CONTINUE 

CALL VSRTR(SPORD(1,1),108,IPERM(1,1)) 
CALL BEIUGR(SPORD(1,I),108,1IOPT,STAT, IER) 
STRES (isl) aSTARC1) 
STRES(2,1)=SOQRT(STAT(5)) 
STRES(3,1)=STAT(4) 

STRES(4,1)=0.0 

STRES(5,1)=SPORD(108,I1) 
STRES(6,1)=SPORD(1,I) 
WRITE(6,70)(STRES(U,1),U=1,6) 

CONTINUE 


ESTIMATE SIGX AND NORMALIZE SPECTRA 


SIGX=(1.+1.3/SQRT(FLOATCINT)))* 
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-SQRT(STRES(4,3)*STRES(4,4)/(2.*INT)) 


DO 230 I1=1,2 

DO 220 J=1, 108 
SPORD(JU,1I)=SPORD(JU,1)/SIGX 
THETA(I )=STRES(3,1)/SIGX 


INITIALIZE DELETION VECTOR 


DO 240 J=1,108 
IDV(J)=1 
ITDEL=O 
NMAX=O 
NDEL(1)=0 
NDEL(2)=0 


FIND NUMBER OF REMAINING CHANNELS AND STD DEV 


NCHAN=O 

SUM(1)=0.0 

SUM(2)=0.0 

DO 260 J=1,108 

IF(IDV(JU).EQ.0)GO TO 260 

NCHAN=NCHAN+ 1 

SUM(1)=SUM(1)+(SPORD(uU, 1)-THETA(1))**2 
SUM(2)=SUM(2)+(SPORD(J,2)-THETA(2))**2 
CONTINUE 

IF (NCHAN.GT.O)GO TO 270 

SDEV(1)=0.0 

SDEV(2)=0.0 

GO TO 275 

SDEV(1)=SQRT(SUM( 1) /NCHAN) 
SDEV(2)=SQRT(SUM(2)/NCHAN) 

DO 280 J=1,4 

ICRIT(U)=0 


CHECK CRITERIA FOR END OF DELETIONS 


IF(SDEV(1).LE.1.25)ICRIT(1)=1 
IF(SDEV(2).LE.1.25)ICRIT(2)=1 
IF(NDEL(1).LT.NMAX )ICRIT(3)=1 

IF (NDEL(2).LT.NMAX)ICRIT(4)=1 
IF((ICRIT(1)+ICRIT(2)).EQ@.2)GO TO 400 
IF((ICRIT(3)+ICRIT(4)).EQ.2)GO TO 400 
ITDEL=ITDEL+1 

NMAX=NMAX+5 

ITER=0O 


REINITIALIZE DELETION VECTOR 


DO 290 J=1,108 
IDV(J)=1 


START DELETIONS WITH SPECTRUM HAVING LARGEST STD DEV 


T=2 
IF (SDEV(2).GT.SDEV(1))I=1 
I=MOD(I,2)+1 


FIND MAXIMUM OUTLIER 
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NDEL(I)=0 

DO 310 J=1, 108 

K=J 

IF(IDV(K).NE.0)GO TO 320 
CONT INUE 

GO TO 370 
SMIN=THETA(I)-SPORD(K,I1) 
DO 330 J=1, 108 

L=109-uJ 
IF(IDV(L).NE.O)GO TO 340 
CONT INUE 

GO TO 370 
SMAX=SPORD(L,1I)-THETA(I) 
IF(SMAX.GE.SMIN)GO TO 350 
SMAX=SMIN 

L=K 


CHECK IF MAX OUTLIER EXCEEDS DELETION CRITERION 


IF(SMAX.LT.6.0)GO TO 370 
NDEL(I)=NDEL(I)+1 


DELETE OUTLIER PLUS ADJACENT CHANNELS 


IDV(L)=0 

IORIG=IPERM(L,I) 

NADU=1 

IF (SMAX.GE.10. )NADJ=2 

DO 360 J=1, 108 
IF(IABS(IPERM(JU,1I)-IORIG).LE.NADJ)IDV(UJ)=0 
CONT INUE 


DELETE MORE IF MAXIMUM ALLOWED NOT REACHED 


IF(NDEL(I).LT.NMAX)GO TO 300 
ITER=ITER+1 


CHECK IF BOTH SPECTRA DONE 
IF(MOD(ITER,2).EQ.1)GO TO 295 
TOLD(1)=THETA(1) 
TOLD(2)=THETA(2) 

FIND HUBER ESTIMATES 


CALL MEST(2) 


CHECK IF MAX NUMBER OF ITERATIONS ALLOWED EXCEEDED 


IF(ITER.GE.6)GO TO 250 

CHECK IF ESTIMATES HAVE CONVERGED 
IF(ABS(THETA(1)-TOLD(1)).GT.0.1)GO TO 285 
IF (ABS(THETA(2)-TOLD(2)).GT.O.1)GO TO 285 
GO. TO 250 

WHEN DELETIONS COMPLETED, WRITE RESULTS 


AMEAN(1)=THETA(1)*SIGX 
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AMEAN(2)=THETA(2)*SIGX 

SDEV(1)=SDEV(1)*SIGX 

SDEV(2)=SDEV(2)*SIGX 

WRITECG 5410) TTDEL, NCHAN, GICRIT(1)41=1.4), SIGX, 


2CAMEAN CI ), SDEV( I) 7 1=1> 2) 


BORMAT CA 2 12 txXe los Ax 40 1.5 (2X, 6102-4 J.) 

FIND BIWEIGHT ESTIMATES 

CALL MEST(1) 

FIND NO. OF OUTLIERS AND TOTAL ABSOLUTE POWER DELETED 


DO 440 I=1,2 

NOUT(I)=0 

SMIN=THETA(I)-5.0 

SMAX=THETA(1I)+5.0 

SUM(I)=0.0 

DO 430 J=1, 108 

IF(IDV(JU).EQ.0)GO TO 420 
IF(SPORD(JU,I).LE.SMIN)GO TO 420 
IF(SPORD(JU,1I).LT.SMAX)GO TO 430 
NOUT(I)=NOUT(I)+1 

SUM(I )=SUM(I )+ABS(SPORD(JU,1)-THETA(I)) 
CONT INUE 

CONT INUE 

SUM(1)=SUM(1)*SIGX 

SUM(2) =SUM(2)*SIGX 
AMEAN(1)=THETA( 1) *SIGX 
AMEAN(2)=THETA(2)*SIGX 

WRITE(6,450) (AMEAN(I),NOUT(I),SUM(I),I=1,2) 
FORMAT (7page2 CE 10842 2X 1352x%) E1024, 2x)) 
GO TO 20 

CONT INUE 

RETURN 

END 


SUBROUTINE TO PERFORM ROBUST MAXIMUM LIKELIHOOD ESTIMATION 
SUBROUTINE MEST(NPSI) 

REAL PSI(600,2),DER(600, 2), THETA(2),SPORD( 108, 2) 
INTEGER IDV( 108) 

COMMON PSI,DER, THETA, SPORD, IDV 
DORZORL=152 

ITER=O0 

T=THETA(I1) 

SUMM=0.0 

SUMD=0.0 

DO 10 J=1,108 

IF(IDV(JU).EQ@.0)GO TO 10 
U=(SPORD(JU,1I)-T)*100. 
INDEX=IFIX(ABS(U)+0.5)+1 

IF (INDEX .GT.600) INDEX=600 
SUMM=SUMM+SIGN(PSI(INDEX,NPSI),U) 
SUMD=SUMD+DER( INDEX ,NPSI) 
CONTINUE 

IF(SUMD.LT.1.0)SUMD=1.0 
DELTA=SUMM/AMAX1(SUMD,5.0) 
T=T+DELTA 

ITER=ITER+1 

IF(ITER.GT.10)GO TO 20 
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361 IF(ABS(DELTA).GT.0.01)GO TO 5 

362 20 THETA(I)=T 

363 30 CONT INUE 

364 RETURN 

365 END 

366 C 

367 @ SUBROUTINE TO FIND DERIVATIVE OF PSI FUNCTION 
368 SUBROUTINE DERIV(PSI,DER) 

369 REAL PSI(600) ,DER(600) 

370 H=0. 1000E-01 

371 DER(1)=(PSI(2)-PSI(1))/H 

372 DER(600)=(PSI(600)-PSI(599))/H 
373 DO 10 I=2,599 

374 10 DER(I )=(PSI(I+1)-PSI(1I-1))/(2.*H) 
375 RETURN 

376 END 


END OF FILE 
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Appendix 5 - FFT Processor Schematics 


The schematics and IC lists are organized in accordance with 
tnegassociatedgcincuit boards. [he circuit] boards and their 

designating letters are the Window Board (W), the Butterfly 

Processor (F), the Power Computation Board (P), and the 


Accumulation Board (A). 
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w-5 


Window Control and Clocking 


Figure A5.3. Schematic W-4 
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Figure A5.5. Schematic F-2 
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Figure A5.10. Schematic F-7 
Input Buffer and FFT Timing 
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Schematic F-8 
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Figure A5.12. Schematic P-2 
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